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1.  FUNDAMENTALS  OF  FLUID  MOTION 

1.1  Introduction 

This  chapter  presents  a  brief  review  of  fluid  flow  fundamentals  pertinent  to  turbulence.  We  expect  them 
to  be  familiar  to  the  reader,  who  may  find  our  particular  viewpoints,  emphasis,  and  compact  notation  helpful 
and  interesting. 

We  will  make  extensive  use  of  the  cartesian  tensor  summation  convention,  where  repeated  indices  imply 
that  the  terms  containing  them  must  be  summed  over  all  possible  coordinate  indices.  An  overdot  ()  will 
be  used  denote  a  partial  derivative  with  respect  to  time,  and  a  subscript  after  a  comma  will  denote  partial 
differentiation  with  respect  to  the  indicated  coordinate  direction;  for  exalnple, 
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“ill  -  “til  +“2i2  +  “3,3  ■ 


We  will  also  use  the  isotropic  tensors  fiy  and  ,  defined  by 


6„  =  (  1  =  i 

1 0  otherwise. 


{1  if  ijk  is  from  the  sequence  123123 
-1  if  ijk  is  from  the  sequence  321321 
0  otherwise 

Various  contractions  will  be  used  frequently,  including 


Sa  —  3  tijkCipq  —  ^jp^kq  ^jq^kp- 

Tensors  are  entities  that,  in  addition  to  being  an  array  of  elements  identified  by  their  subscripts,  trans¬ 
form  in  a  very  special  way  when  the  coordinate  system  is  transformed  by  rotation.  A  tensor  that  is  totally 
unchanged  by  an  arbitrary  rotation  of  the  coordinate  system  is  called  isotropic.  Any  second-order  isotropic 
tensor  must  be  a  scalar  times  S, y,  and  any  third-order  isotropic  tensor  must  be  a  scalar  times  f,y*.  Moreover, 
any  higher-order  isotropic  tensor  must  be  expressible  in  terms  of  the  various  possible  combinations  of  these 
two  tensors,  and  hence  they  are  fundamental  building  blocks  in  all  sorts  of  physical  modeling,  including 
viscous  flow  and  turbulence. 

1.2  The  basic  equations 

The  basic  equations  are  derived  by  application  of  basic  principles  to  an  elemental  control  volume  (Fig. 
1.2.1).  The  conservation  of  mass  gives 

P  +  (e“j).j  =  0  (1.2.1) 

where  p  is  the  fluid  density,  and  uy  is  the  fluid  velocity  component  in  the  direction.  This  is  also  called 
the  continuity  equation. 
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Figure  1.2.1  Control  volume  for  basic  equation  derivation 

The  momentum  equation  is 

(pu,)  +  (ptiytiy),y  =  <Ty,,y  +/,  (1.2.2) 

where  Oji  is  the  stress  in  the  »,k  direction  on  a  control  volume  surface  perpendicular  to  the  jth  axis,  and  /, 
is  the  body  force  (per  unit  volume)  in  the  jth  direction. 

The  conservation  oj  energy  requires  that 


(P‘ o)  +  (pujeo),,  =  (<ry;u,)»  +/.U,  -  gy.y.  (1.2.3) 

Here  eo  =  e  +  ^ V 3  is  the  total  energy  per  unit  mass,  where  e  is  the  internal  energy  per  unit  mass  and 
V1  —  u, u, ,  and  </y  is  the  conduction  heat  flux  (flow  rate  per  unit  area)  in  the  j,h  direction  outward  from  the 
elemental  control  volume.  The  first  term  on  the  right  represents  the  power  input  by  the  surface  forces  per 
unit  volume,  and  the  second  that  by  the  body  forces. 

The  entropy  balance  is 

(ps)  +  (puy s),y  =  q>  -  (?y/T),y  (1.2.4) 

where  s  is  the  entropy  per  unit  mass,  T  is  the  absolute  temperature,  and  <p  is  t.ie  entropy  production  rate 
per  unit  volume.  Here  the  term  gy/T  represents  the  entropy  flux  associated  with  the  heat  flux  gy.  The 
second  law  of  thermodynamics  requires  that  the  entropy  production  be  non-negative, 


ip  >  0. 


(1.2.5) 


These  ideas  are  useful  in  assessing  constitutive  models  for  the  stress  tensor  and  heat  flux  vector,  and  in 
identifying  the  processes  that  produce  entropy  (dfaiipole  energy)  in  viscous  flows. 

1.3  The  stress  tensor 


The  stress  tensor  <r,y  must  be  symmetric.  This  fact  can  be  established  by  performing  a  moment  of 
momentum  analysis  on  the  elemental  control  volume  of  Fig.  1.3.1.  The  torques  of  the  stress  terms  are  all 
of  order  dijdij,  and  the  moments  of  the  momentum  flows  and  body  forces  are  all  of  higher  order,  hence 

Ou  =  <7ji. 
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Figure  1.3.1  Control  volume  for  stress  tensor  symmetry  derivation 


The  tensor  can  be  split  into  two  parts: 


<7,y  =  -PSij  +  T.y.  (1.3.1) 

The  P  term  represents  the  isotropic  component  of  the  (inward)  normal  stress;  r; y  is  the  deviations  from  this 
isotropic  stress,  attributed  phenomenologically  to  viscosity.  From  a  molecular  point  of  view,  cr,y  arises  from 
molecular  transport  of  momentum;  the  isotropic  part  P  is  determined  by  the  average  using  the  probability 
distribution  for  molecular  velocities  (e.g.  Boltimann),  and  r,y  arises  from  anisotropy  in  the  probability 
distribution. 
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1.4  Thermodynamic  properties  and  concept* 

The  internal  energy  e  reflect!  the  randomly-oriented  energy  of  molecular  translation,  rotation,  vibration, 
and  other  microscopic  energy  modes  (chemical  bonds,  etc.).  In  general,  e  depends  upon  the  thermodynamic 
state  (i.e.  T  and  p),  but  for  idealised  gas  and  incompressible  liquid  models  depends  only  on  temperature. 
It  is  sometimes  called  the  thermal  energy,  and  is  all  too  frequently  confused  with  heat  (qy),  which  is  the 
transport  of  energy  by  disordered  molecular  processes.  The  internal  energy  of  an  object  can  be  increased  by 
a  transfer  of  energy  as  either  heat  or  work,  and  the  energy  flowing  as  heat  can  come  either  from  a  source  of 
internal  energy  01  mechanical  energy  (kinetic,  potential,  or  work).  The  internal  energy  is  a  thermodynamic 
property  of  matter,  the  heat  transfer  is  not.  The  confusion  between  heat  and  internal  energy  is  an  infortunate 
remnant  of  the  caloric  theory  of  heat,  but  perhaps  understandable  since  the  theory  was  discarded  only  about 
a  century  ago. 

The  entropy  can  be  thought  of  as  a  measure  of  the  degree  of  randomness  at  the  molecular  level,  and 
in  modern  thermodynamic  treatments  the  tempet  ature  is  interpreted  as  a  measure  of  the  sensitivity  of  this 
randomness  to  changes  in  energy  at  constant  density.  Orderly  microscopic  exchanges  of  energy  (e.g.  as  work 
or  as  bulk  kinetic  energy)  have  no  associated  entropy  transport.  But  heat,  the  microscopically  disordered 
transport  of  energy,  does  carry  entropy  with  it,  and  it  may  be  shown  that  this  entropy  transfer  flux  is  qy/T. 
For  more  discussion  of  these  important  thermodynamic  concepts  from  this  viewpoint,  see  Reynolds  and 
Perkins  (1977). 

It  is  usually  assumed  that  as  far  as  the  thermodynamic  properties  are  concerned  the  fluid  is  in  a  state 
of  local  equilibrium,  and  hence  that  the  usual  relations  between  thermodynamic  properties  are  valid.  Thus, 
the  Gibbs  equation  is  used  to  relate  entropy  changes  to  energy  and  density  changes, 


Tds  =  ie  +  Pd(l/p).  (1.4.1) 

The  enthalpy  h  is  defined  as 

h  =  e  +  P/p  (1.4.2) 

and  represents  the  sum  of  the  convected  internal  energy  and  flow  work  associated  with  the  transport  of 
a  unit  mass  of  fluid  across  a  control  volume  boundary.  We  emphasize  that  it  is  the  internal  energy  that 
appears  in  the  basic  energy  balance  equation. 

An  alternative  form  of  the  energy  equation  is  obtained  using  (1.3.1)  in  (1.2.3),  moving  the  pressure  term 
to  the  left  hand  side; 

(p'o)  +  (pUjho),j  =  (rjiU,),y  +/,u,  -  qy,y.  (1.4.3) 

Here  h0  -  h  +  1 V 2  is  the  stagnation  enthalpy.  Note  that  the  enthalpy  appears  as  the  convected  energy  per 
unit  mass  (internal  energy  e  plus  flow  work  P/p),  but  the  internal  energy  e  appears  in  the  energy  storage 
rate  term.  A  common  error  is  the  use  of  enthalpy  in  both  places. 

1.6  Kinematics  of  motion 

Any  deformation  rate  u;,y  can  be  decomposed  into  the  sum  of  a  strain  rate  S;y  and  a  rotation  rate  fl;y, 


u<>/  —  2  ’>  )  +  2  (“*’*  ui”  1  —  +  fi,y. 


(1.5.1) 


Note  that  the  strain  rate  is  a  symmetric  tensor  and  the  rotation  rate  is  antisymmetric.  They  play  quite 
different  roles  in  fluid  mechanics,  particularly  in  turbulence,  and  for  this  reason  we  prefer  forms  of  the 
equations  that  make  their  presence  or  absence  very  clear. 

1.6  Mechanical  and  thermal  energy  equations 

The  fundamental  equations  may  be  combined  to  derive  an  equation  describing  the  transport  of  macro¬ 
scopic  mechanical  energy  and  another  describing  the  transport  of  internal  energy.  The  mechanical  energy 
equation  is  derived  by  contracting  the  momentum  equation  with  the  velocity;  multiplying  (1.2.2);  by  Uj, 


(p\V*)  +  (pu]\v7)'J  =  +  At',. 


The  right  hand  side  may  be  written  as 


(1.6.1) 


*1 


)  ,y  (TijUi ,y  +  tig . 


Then,  using  (1.3.1)  to  split  the  stress  tensor  and  (1.5.1)  in  place  of  velocity  gradients,  and  noting  that 
r,yfl,y  =  0  since  r,y  is  symmetric  and  fi.y  is  antisymmetric  (sum  over  both  repeated  indices  is  implied), 
(1.6.1)  can  be  rewritten  as 

(/)-Va)  +  (puy -V3)  ,y  =  PSjj  —  (Pu,),y  +  /.u,  +  (r,yu,),y  -  (r,y£;y).  (1.6.2) 

The  sum  on  the  right  represents  the  input  of  macroscopic  mechanical  energy  to  the  control  volume, 
which  shows  up  as  an  increase  in  the  kinetic  energy  of  the  flow.  Two  of  these  terms  appear  as  power  inputs 
in  the  thermal  energy  equation  (1.6.4)  but  with  opposite  sign,  and  hence  these  terms  represent  exchanges 
between  thermal  and  mechanical  energy.  The  first  is  PSjj,  which  is  the  rate  of  energy  transfer,  per  unit 
volume,  from  thermal  energy  to  mechanical  energy  due  to  expansion  of  the  fluid.  The  second  is  which 

represents  the  transfer  of  mechanical  energy  to  thermal  energy  by  viscous  forces.  This  is  the  only  viscous 
term  involved  in  the  entropy  production  equation  (1.7.3),  and  hence  this  is  the  only  viscous  term  properly 
termed  diiipation.  Since 

I  (rlyu.),yd3x  =  0  (1.6.3) 

if  the  integral  is  taken  over  a  volume  where  either  the  velocity  or  stress  is  zero  on  the  boundaries,  this  viscous 
term  has  no  global  effect;  it  represents  reversible  viscous  power  input  to  the  control  volume  from  surrounding 
fluid.  The  term  containing  ft  is  the  power  input  from  body  forces,  and  the  (P Uj),y  term  represents  power 
output  by  flow  work. 

The  thermal  energy  equation  is  obtained  by  subtracting  the  mechanical  energy  equation  (1.6.2)  from 
the  total  energy  equation  (1.2.3),  and  is 

(p<)  +  (ptiye).y  =  -PSjj  +  (r.yS.y)  -  gy.y  .  (1.6.4) 

Here  PSjj  represents  the  power  output  from  thermal  energy  due  to  expansion  of  the  fluid,  TtJ S„  is  the  power 
input  to  the  thermal  energy  due  to  irreversible  viscous  effects,  and  gy„  is  the  net  power  output  due  to  heat 
conduction,  all  per  unit  volume. 

Note  that  the  enthalpy,  which  appeared  in  the  alternate  form  of  the  total  energy  equation  (1.4.3),  does 
not  appear  in  the  thermal  energy  equation.  We  have  derived  the  thermal  energy  equation  correctly  from 
first  principles.  One  must  be  wary  in  reading  literature  where  the  thermal  energy  equation  is  developed  from 
a  “heat  balance’,  because  there  is  no  such  principle  as  the  conservation  of  heat. 

1.7  Irreversibility  rate  equation 

Using  the  conservation  of  mass  equation  (1.2.1)  in  the  Gibbs  equation  (1.4.1), 

pTDs  =  pDe  +  PSa  (1.7.1) 

where  D  denotes  the  substantial  derivative 

£>()  =  ()  +  Uj(),j.  (1.7.2) 

Using  the  thermal  energy  equation  (1.6.4)  and  the  entropy  balance  (1.2.4),  this  yields  an  expression  for  the 
irreversibility  rate, 

Ty>=  -jiqjT,j+rijS,j>0.  (1.7.3) 

This  clearly  identifies  the  viscous  dissipation  term  as  discussed  above,  and  provides  a  neat  framework  for 
evaluation  of  consitutive  models  for  the  heat  flux  or  viscous  stresses. 

1.8  Constitutive  equations 

The  theory  of  linear  algebra  is  extremely  helpful  in  developing  constitutive  models  for  the  heat  flux  and 
viscous  stresses,  and  also  for  developing  turbulence  models.  We  will  use  these  ideas  to  review  the  constitutive 
equations  so  as  to  set  the  stage  for  later  use  of  these  ideas  in  developing  turbulence  models. 

The  most  general  vector  /,  that  is  a  function  of  only  one  other  vector  ti<  is 

fi  =  Cvi  (18.1) 

where  the  coefficient  C  can  be  a  function  of  scalars,  including  the  invariant  of  the  vector  (its  magnitude 
VfcVfc).  Higher-order  terms,  such  as  u;u*vj„  need  not  be  added  since  they  are  represented  by  allowing  the 
coefficient  to  depend  on  the  invariant  of  v.  Thus,  if  one  assumes  that  the  heat  flux  vector  q,  is  a  function  of 
the  temperature  gradient  vector  T„,  the  most  general  form  is  the  familiar  Fourier  heat  conduction  law, 
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<b  =  ~kT,i  (1.8.2) 

where  k  is  the  tnermal  conductivity ,  which  depends  to  first  approximation  on  the  temperature  and  may, 
in  mgher  approximations,  also  depend  on  the  scalar  TlkT,k.  It  is  generally  believed  that  (1.8.2)  describes 
heat  conduction  in  fluids,  except  perhaps  in  regions  of  very  strong  temperature  gradient  such  shock  waves 
or  combustion  fronts. 

The  most  general  tensor  at)  that  is  a  function  of  only  one  other  tensor  b, ,-y  is,  in  three  dimensions, 

o;y  =  ASij  +  Bbij  +  Cbjj  (1.8  3) 

wl^re  —  bikbkj.  The  coefficients  A,  B  and  C  may  depend  on  relevant  scalars,  including  the  three  scalar 
invariants  of  th'  tensor  b.  Higher-order  terms  such  as  i*y  =  Jj* 6ty  need  not  be  added  since,  by  the  Cayley- 
Ilamilton  thicrem,  they  can  be  expressed  in  terms  of  lower-order  terms  and  the  invariants  of  b  and  hence 
are  already  included  in  (1,8.3).  Therefore,  if  it  is  assumed  that  the  viscous  stress  tensor  r^y  is  a  function  of 
the  local  strain-rate  tensor  5,y,  this  functional  dependence  must  be  of  the  form 

r.y  =  A6„  +  BS„  +  CS?  (1.8.4) 

where  the  coefficients  may  depend  on  scalars,  such  as  temperature,  density,  or  the  invariants  of  S.  This  is 
called  the  Stokes  model  for  viscous  stresses. 

The  rms  strain-rate  S  —  (SiySy,)1^3  is  a  reciprocal  time  scale  for  the  fluid  deformation.  If  this  time 
is  long  compared  to  molecular  collision  times,  then  the  strain  is  considered  weak  and  only  linear  terms  in 
(1.8.4)  are  used.  This  leads  to 

rtJ  =  My  +  BSi,  (1.8.5) 

where  A  can  depend  at  most  linearly  on  the  invariants  of  S,  and  B  must  be  independent  of  S.  If  it  further 
assumed  that  P  =  -  Jo,-,,  then  by  (1.3.1) 


Tkk  =  0  =  3.4  +  BSkk 


so 

A  =  -^BSkk.  (1.8.6) 

for  a  simple  sheai  ing  flow  where  the  only  non-tero  strain-rate  elements  are 


c  _  -  1  dui 

one  defines  the  fluid  viscosity  ft  by 

Tyi  =  2ftSn  (1.8.7) 

from  which,  using  (1.8.5),  it  follows  that  B  =  2 ft.  The  resulting  Newtonian  constitutive  equation  is 


2 

Uj  —  2ftSij  ~  ~ftSkk6ij.  (1.8.8) 

Note  that  the  Newtonian  constitutive  equation  assumes  only  that  the  viscous  stress  tensor  is  a  trace-free 
linear  function  of  the  local  strain  rate;  this  assumption  is  believed  to  be  quite  adequate  for  many  continuum 
fluid  flows.  The  modal  fails  in  strong  shock  waves  (normal  stresses  are  incorrect)  and  in  flow  of  polymers 
(rotation  rates  are  also  important). 

Using  (1.8.2)  and  (1.8.8)  in  (1.7.3),  the  irreversibility  rate  becomes 

3V  =  fTti  Tti  +2yj(5,y5,y  -  ^S„Skk)  >  0,  (1.8.9) 

It  is  clear  that  the  heat  flux  term  is  positive-definite.  It  is  left  as  an  exercise  to  demonstrate  that  the 
strain-rate  term  is  also  positive-definite  (Hint:  evaluate  in  the  principal  coordinates  of  Sis  by  expressing  the 
diagonal  elements  as  the  components  of  a  vector  in  polar  coordinates). 

1.9  Vortlclty 

Vorticity  is  one  of  the  most  fundamental  concepts  in  fluid  mechanics,  and  probably  the  most  important 
concept  in  turbulence.  The  vorticity  vector  u.  is  defined  by 


<0,  =  C,y*Ufc,y  . 


(1.9.1) 
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Note  that  =  0  and  hence  the  vorticity,  by  definition,  ii  divergence  free, 

Ui»=  0. 

By  the  definition!,  the  vorticity  ii  related  to  the  rotation  rate, 

w,  = 

Taking  ep,,'  x (1,9. 3)<  ,  one  find! 

1 

fl,p  —  2epi‘u>- 


(1.9.2) 

(1.9.3) 

(1.9.4) 


finds 


Th«  vorticity  field  can  be  thought  of  as  contributing  to  the  velocity  field.  Forming  fPqiX(l.9.l)ilfll  one 


—  ^pqi^ijkuk$jq  —  uqipq 


Vifkk  —  tikjWjik  "b(ufc»fc  )»i  • 


(1.9.5) 


Thii  ii  a  Poisson  equation  for  the  velocity,  analogue  to  the  equation  for  temperature  in  a  heat  conducting 
medium  with  distributed  eourcee.  Eqn.  (1.9.5)  displays  two  ‘sources*  of  velocity,  namely  vorticity  (more 
specifically  vorticity  gradients )  and  flow  divergence  (expansion  or  compression).  In  addition  to  the  velocity 
generated  by  these  sources,  one  can  also  have  an  additional  component  of  velocity  satisfying  the  Laplace 
equation  u ,-,**=  0.  FYom  (1.9.5)  we  see  that  this  component  could  be  thought  of  as  arising  from  uniform 
vorticity  (a  solid-body  rotation)  and  uniform  irrotational  expansion,  of  which  irrotational  flow  at  constant 
density  is  a  special  case. 

The  part  of  the  velocity  field  due  to  the  vorticity  gradients  may  be  found  using  the  general  solution  to 
the  Poisson  equation;  at  any  instant  in  time,  this  solution  is 


“t(x)  =  ~  J  G(x,x')(ik, u>j,k  (x')dV 


(1.9.6) 


where  (.(x.x')  is  the  Green’s  /unction  for  the  Poisson  solution  in  the  flow  domain,  and  d3x'  represents  an 
element  of  volume  for  the  interation  over  the  flow  domain.  The  Green’s  function  for  an  infinite  domain  is 


C(x.x')  = 


-1 


4»n/(*»  -  <)(*«  -  *;) 


(1.9.7) 


Using  this  Green’s  function  in  (1.9.6),  and  integrating  by  parts  to  transfer  the  k  differentiation  from  the 
vorticity  to  the  Green’s  function,  one  finds 


«.(*)  =  j 


‘ ilk 


[*k  ~  x'k) 


4*|(zn  -  z'J(x„  -  z')p/’ 


Uj<Px‘ 


(1.9.8) 


This  is  called  the  Biot-Savart  equation.  It  gives  that  portion  of  the  velocity  field  arising  from  vorticity,  for  an 
infinite  flow  domain.  Computational  methods  in  which  markers  track  the  motion  of  vorticity-bearing  fluid 
use  the  Biot-Savart  equation  to  compute  the  velocity  field;  this  is  an  efficient  calculation  if  the  vorticity  is 
highly  concentrated  and  most  of  the  fluid  has  negligible  vorticity,  and  there  are  many  interesting  problems 
in  turbulence  that  can  be  addressed  in  this  manner. 

We  emphasize  that  all  of  the  features  of  vorticity  discussed  thus  far  are  kinematic  in  nature,  and  apply  in 
either  compressible  or  incompressible  flows.  In  the  next  section  we  will  adress  the  dynamics  of  the  vorticity. 

1.10  Vorticity  dynamics 

Using  the  continuity  equation  (1  2.1)  and  the  stress  tensor  split  (1.3.1),  the  momentum  equation  (1.2.2) 
can  be  written  as 


Vjk  "f"  UqU/tiq  —  J*ik  "f*  fk]  • 


(1.10.1) 


Taking  (1.10.1)  k,j  one  obtains 


+  Uy Wj  ,y  —  [“  {^kqyq  P| fc  “h A  )  ] , 


Using  (1.5.1)  and  (1.9.4),  the  first  term  on  the  right  is  exactly 

WqSiq  —  SqqU>i  • 
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The  pressure  term  can  be  expanded  into  two  terms,  one  of  which  vanishes,  and  the  vorticity  equation  becomes 

o!i  +  =  UjSij  -  SjjUii  +  tijk  [(^7t,„),j  +^P,k  p, j  +(~)u]  .  (1.10.2) 

Note  that  the  left-hand  side  of  (1.10.2)  represents  the  rate  of  change  of  vorticity  following  a  fluid  particle. 
Thus,  the  terms  on  the  right  display  the  processes  that  can  give  rise  to  changes  in  vorticity  of  a  fluid  particle. 
The  first  term  represents  the  straining  of  vortex  filaments,  and  is  a  crucial  term  in  turbulence;  in  a  two- 
dimensional  flow,  this  strain  is  always  in  planes  perpendicular  to  the  vorticity,  and  hence  there  is  no  vortex 
stretching  in  two-dimensional  flow.  The  second  term  shows  that  fluid  compression  (S**  <  0)  tends  to  amplify 
the  voiticity,  and  expansion  to  attenuate  it.  The  term  containing  q,,  represents  viscous  effects,  including 
diffusion.  The  term  containing  pressure  gradients  and  density  gradients  shows  that  these  may  combine  to 
act  as  a  source  for  vorticity,  if  these  gradient  vectors  have  a  non-sero  cross  product;  this  term  is  important 
in  the  atmosphere.  Body  force  gradients  can  also  generate  vorticity;  but  body  forces  are  often  conservative, 
i.e.  of  the  form 

fk  =  p4‘,k  (1.10.3) 

where  <J>  is  a  scalar  potential,  and  (1.10.2)  shows  that  such  forces  do  not  generate  vorticity. 

In  a  Newtonian  flow  where  p  =  p(t),  n  is  constant,  and  /*  =  p<t>,k,  the  vorticity  equation  becomes 


Ui  +  UjUi,j  =  uijSi,-  -  UiSjj  +  vwujj .  (1.10.4) 

This  is  the  form  to  which  we  will  refer  most  often  in  our  studies  of  turbulence;  it  emphasizes  the  interaction 
between  strain-rate  and  vorticity  that  is  so  important  in  turbulence. 

One  usually  sees  the  vorticity  equation  with  the  first  term  on  the  right  in  (1.10.4)  replaced  using  an 
identity  derived  from  (1.5.1)  and  (1.9.4), 

WyS,,  =  w,u,„,  (1.10.5) 

We  prefer  ( 1.10.4)  because  it  makes  the  interaction  between  vorticity  and  strain-rate  very  clear. 

1.11  Vortex  line*  and  tubes 


/ 
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Figure  1.11.1  Velocities  along  a  vortex  line 
A  vortex  line  is  a  line  everywhere  tangent 


to  the 

Sxi_ 

St 


W.  vector.  Along  a  vortex  line  (see  Figure  1.11.1) 


W 


(1.11.1) 


Vortex  lines  move  as  the  fluid  moves.  For  inviscid,  incompressible  flow,  the  vortex  lines  move  with  the 
fluid.  This  fact  is  extremely  helpful  in  understanding  fluid  flows  in  general  and  turbulence  in  particular,  and 
forms  the  basis  for  an  important  class  of  numerical  methods  for  simulating  turbulent  flow. 

We  will  now  prove  this  important  fact  about  vortex  lines.  Let  w  be  the  vorticity  at  the  center  of  an 
elemental  segment  Ss  of  a  line  marked  in  the  fluid  along  a  vortex  line  at  time  t.  The  rate  of  change  of  the 
vorticity  following  the  fluid  particle  attached  to  the  center  of  this  line  given  by  (1.10.4).  Neglecting  the 
viscous  term,  and  assuming  constant  density  (so  that  5**  =  0),  and  using  (1.10.5),  the  rate  of  change  of  the 
vorticity  of  this  fluid  particle  is 

a/,=t5,u,„.  (1.11.2) 

The  right  hand  side  of  (1.11.2)  is  evaluated  at  time  t  using  (1.11.1)  to  express  Sii  in  terms  of  Ss,  yielding 


(1.11.3) 
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Since  each  end  of  the  line  moves  with  iti  own  velocity, 

(1.11. 

We  now  examine  the  changes  in  6x i/u»i : 

d  .in,  1  d{6xi)  Sxi  dZj  ,)u 

dt '  Ui  tii  dt  u7  dt 

Uiing  (1.11.4)  for  the  first  term  on  the  right  and  (1.11.3)  for  the  second,  the  right  hand  side  is  tero,  and 
hence  Sxi/iiii  is  constant  in  time.  The  same  is  true  for  the  other  two  components.  Hence,  if  the  line  was 
originally  a  vortex  line,  it  will  remain  a  vortex  line,  as  we  have  claimed.  Moreover,  since  Sx i  =  CjyC,  it 
follows  that  |{*i|  will  be  proportional  to  the  line  length. 

We  can  form  a  vortex  tube  from  a  set  of  nearby  vortex  linos.  In  inviscid,  incompressible  flow,  this 
tube  will  move  with  the  fluid,  and  can  be  stretched  by  strain  along  its  length.  This  strain  will  intensify  the 
vorticity  in  the  tube.  Since  the  fluid  is  incompressible,  and  the  tube  is  imbedded  in  the  fluid,  stretching  the 
tube  reduces  its  diameter.  The  increase  in  vorticity  can  be  though  of  in  terms  of  the  increased  rotational 
rate  necessary  to  maintain  conserved  angular  momentum  as  the  tube  decreases  in  diameter.  These  processes 
of  vortex  convection  and  stretching  by  the  flow  are  central  in  turbulence. 

It  is  left  as  an  exericise  to  show  that,  in  inviscid,  compressible  flow,  lines  everywhere  tangent  to  &/p 
move  with  the  fluid.  This  fact  may  be  useful  in  simulations  of  compressible  turbulence. 

2.  TURBULENCE  EQUATIONS 


d(Sx<) 

dt 


2.1  Averaging  concepts 

Different  kinds  of  averaging  procedures  are  appropriate  for  different  situations.  In  situations  that  are 
statistically  steady,  the  time  average  is  useful.  Denoting  a  random  field  by  /(x,t),  its  time  average  is 

7(x)  =  lim  i  f  /(x,  t)dt.  (2.1.1) 

T-.00  l  Jo 

The  time  average  can  not  be  used  in  fields  that  are  statistically  developing  in  time.  But  if  the  field  is 
statistically  homogeneous,  i.e.  statistically  the  same  at  all  space  points,  then  a  volume  average  can  be  used, 

(»■>■») 

If  the  field  is  not  statistically  steady  or  homogeneous  in  space,  but  is  homogeneous  on  planes  or  along 
lines,  averages  on  the  planes  or  lines  can  be  used.  But  if  the  field  is  not  statistically  the  same  in  time  or  any 
space  dimension,  one  has  to  resort  to  the  concept  of  ensemble  averaging,  i.e.  averaging  over  a  large  set  of 
(usually  hypothetical)  similiar  experiments: 


1  N 

/(*.«)  =  Jj  5Z  /nfc'1)- 


(2.1.3) 


One  must  always  be  careful  to  choose  an  averaging  concept  appropriate  to  the  situation.  It  will  be 
assumed  that  an  averaging  process  has  been  chosen  that  commutes  with  differentiations  with  respect  to 
both  time  and  space;  the  ensemble  average  always  has  this  property. 

2.2  Turbulence  decomposition 

Each  variable  in  a  random  field  is  represented  as  the  sum  of  its  average  and  its  fluctuation, 


It  follows  that 


/  =  /+/'. 

(2.2.1) 

tch  that 

7  =  0. 

(2.2.2) 

7h  =  Jh 

(2.2.3) 

Jh'  =  0. 

(2.2.4) 

and 
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In  compressible  flows,  mass-weighted  averaging  is  often  employed.  The  methods  for  doing  this  averaging  are 
simple  extensions  of  those  given  above. 

Most  turbulence  literature  concerns  incompressible  Sows.  However,  there  is  a  class  of  compressible 
flows  that  can  be  handled  as  a  very  modest  extension  of  incompressible  flows,  namely  flows  where  p  =  p(t) 
(uniform  density  flows).  Many  practical  flows  fall  in  this  class,  including  flow  in  an  internal  combustion 
engine  cylinder.  The  equations  for  uniform  density  flow  are  much  simpler  than  those  of  full  compressible 
flow,  and  so  in  the  interest  of  simplicity  much  of  what  follows  will  be  limited  to  uniform  density  flows  with 
constant  viscosity  pi. 

2.3  Governing  equations 

If  p  —  p(t),  the  continuity  equation  reduces  to 

p  +  pui,i  =  0.  (2.3.1) 

We  will  write  the  turbulence  decomposition  with  capital  letters  for  mean  quantities  and  lower  case  letters 
for  the  fluctuations, 

p  =  P- f  p'  (2.3.2a) 

Ui*=Ui  +  uj.  (2.3.26) 

Inserting  these  decompositions  into  the  continuity  equation  (2.3.1),  and  averaging,  we  obtain  the  mean 
continuity  equation 

p+pUiti  =  0.  (2.3.3) 

Subtracting  this  result  from  (2.3.1)  we  obtain  the  fluctuating  continuity  equation, 


=  o 


(2.3.4). 


Note  that,  for  uniform  denity  flow,  the  fluctuation  velocity  Held  is  divergence-free,  as  would  be  the  mean 
velocity  field  if  the  flow  were  incompressible. 

For  uniform  density  flow  with  p  ^constant  and  fk  =  0,  the  momentum  equation  (1.2.2)  reduces  to 

t*i  +  u,u i,j  =  - -p„  +!/«,,„• .  (2.3.5) 

P 

Introducing  the  turbulence  decomposition,  averaging,  and  making  use  of  (2.3.4),  the  mean  momentum  equa- 
t  'on  is  found  as 

Ui  +  UpUi.i  =  - -P,i  +uUi,jj  (2.3.6) 

P 

where 

Rii  =  (2.3.7) 

The  quantity  -pR\j  appears  in  (2.3.6)  like  a  stress,  and  so  is  called  the  Reynolds  stress  after  O.  Reynolds, 
who  intp  luced  the  basic  decomposition. 

Equations  (2.3.3)  and  (2.3.6)  would  permit  calculation  of  the  mean  density  and  velocity  field  if  the 
Reynolds  stresses  were  all  known.  Since  they  are  not  known,  we  have  a  closure  problem,  which  can  be 
addressed,  but  not  solved,  by  further  development  of  the  equations  for  the  Reynolds  stresses. 

An  alternative  way  of  thinking  about  the  turbulence  ‘forces*  has  some  physical  appeal.  From  (1.9.1)  it 
follows  that 


“X  =  “p 


Multiplying  by  e,,p 


c,ipw'u(,  -  u[>  —  (<5pj6?k  6pk6u)ukij  up  -  ui>p  UJ>  u!»«  up- 


Using  (2.3.4),  this  produces 


or  equivalently 


fe.>uiup  =  K“p ).p  -~Kup)-« 


(«!«').  J  =  5K“y)»  +fi/fcwyu*- 


p'-p+rf 


We  define 


(2.3.10) 


MO 


where  the  mean  iquare  fluctuation  velocity  is 


93  =u'uj. 


It  is  convenient  to  denote  the  mean-convection  substantial  derivative  by 

E()  =  ()  +  lM).y 

Then,  using  (2.3.9)  in  (2.3.6),  the  mean  momentum  equation  becomes 


(2.3.11) 


(2.3.12) 


DUj  = 


(2.3.13) 


In  this  alternative  representation  the  turbulence  provides  a  contribution  to  the  normal  stress  in  P’  and  a 
turbulent  body  force,  but  no  shear  stress.  The  potential  of  this  alternative  view  of  turbulence  forces  remains 
to  be  investigated. 

2.4.1  Mean  vorticily  equation 

The  turbulence  decomposition  of  the  vorticity  is 

u)i=ni+ar;.  (2.4.1) 

The  mean  strain  rate  and  rotation  rate  are  denoted  by  S,i  and  Cl{j,  respectively,  and  the  fluctuation  strain 
rate  by  s(y.  By  continuity  for  uniform  density  flow,  s',  =  0.  Introducing  the  turbulence  decomposition  into 
(1.10.4),  and  averaging,  the  mean  vorticity  equation  is  found  to  be 

DOi  =  (IjSi,  -  U,S„  +  I/O,,,,  -K^),3  (2.4.2) 

Note  the  appearance  of  the  turbulence  body  force  term  w'Uy  in  the  equation. 

2.6  Turbulent  fluctuation  equations 

The  fluctuating  continuity  equation  is  (2.3.4).  Subtraction  of  the  mean  momentum  equation  (2.3.6) 
from  the  full  equation  (2.3.5)  gives  the  fluctuating  momentum  equation, 

Du'i  =  -u'jUiti  -(u'u'  -  u'u'),y  -  Ip'„  +i/u',„  .  (2.5.1) 

Subtraction  of  (2.4.2)  from  (1.10.4)  gives  the  fluctuating  vorticity  equation 

Du'i  =  +w'  S(J  -  u>'5„  + 

-u' n.,y  -(u'w'  -  u'w')„  +(w'j'y  -  <j's;y)  +  i/w',„  .  (2.5.2) 

By  takmg  (2.5.1),,-  one  obtains  ar  equation  for  the  fluctuating  pressure, 

~P  lit  =  Uj,i  —  (u,-,y  Uy —  u(.,y  Uy,,')  (2.5.3) 

These  equations  are  useful  for  deriving  equations  relating  statistical  properties  of  the  turbulence  and  in  the 
study  of  the  dynamics  of  turbulent  fluctuations. 

2.6  Kinetic  energy  equations 

The  transport  equation  for  the  turbulent  kinetic  energy 

(2.6.1) 

is  derived  by  multiplying  (2.5.1),  by  u'  and  averaging.  After  some  rearranging,  one  obtains 

=  (2.6.2) 


P  =  -«'u'50 


Ml 


U  the  rate  of  turbulent  energy  production, 


0  = ,y 


(2.6.4) 


if  the  (homogeneous)  rate  of  turbulent  energy  dissipation,  and 


j}.  =  V“' + 


(2.6.5) 


if  the  turbulent  kinetic  energy  flux  in  the  direction. 

Note  that  P  involves  the  product  of  turbulent  stresses  and  mean  strain  rates,  and  that  the  mean 
rotation  rate  does  not  appear  explicitly  in  the  turbulent  kinetic  energy  equation  (though  it  may  influence 
the  turbulent  kinetic  energy  by  altering  terms  in  the  equation).  P  arises  from  the  stretching  of  the  tangle 
of  vortex  filaments  that  make  ip  the  turbulence  by  the  mean  deformation.  P  is  almost  always  found  to  be 
positive,  although  there  are  situations  in  which  it  is  negative. 

Since  the  source  of  turbulent  kinetic  energy  is  the  mean  flow,  the  production  term  should  appear 
with  opposite  sign  in  the  evolution  equation  of  the  mean  kinetic  energy.  Multiplying  (2.3.6),  by  C/, ,  and 
rearranging,  the  mean  kinetic  energy  equation  is 

D(\UiUi)  =  -(-Pf/,)„+-/J5„  -  2i +  2v(UiSij),j -P  -  (2.6.6) 

2  p  P 

Indeed,  —  P  does  appear  on  the  right,  indicating  a  drain  from  the  mean  kinetic  energy.  The  two  pressure 
terms  represent  the  power  associated  with  flow  work  and  the  power  transfer  from  internal  energy  due  to 
expansion  of  the  flow.  The  first  viscous  term  is  the  rate  of  dissipation  of  mean  kinetic  energy  by  viscous 
effects  (see  1.8.9),  and  the  second  is  the  reversible  viscous  power  transfer.  The  last  term,  which  integrates 
to  sero  over  a  large  volume  of  flow  bounded  by  non-turbulent  fluid,  represents  an  internally  conservative 
transfer  of  mean  kinetic  energy  within  this  volume. 

We  have  been  careful  to  call  P  the  homogeneous  dissipation,  because  (as  shown  in  the  next  chapter) 
only  in  homogeneous  turbulence  does  it  represent  the  true  rate  of  energy  dissipation.  FVom  (1.8.9)  the  true 

dissipation  rate  is  _  _  _ 

e  =  I'r'yfJj  =  2«/»|ys{y  =  */u',y  (u',y  ).  (2.6.7) 

The  right  hand  side  of  the  turbulent  kinetic  energy  equation  can  be  modified  to  replace  P  by  e,  with  a 
concurrent  modification  in  the  definition  of  the  flux.  This  is  left  as  an  exercise. 

2.7  Reynolds  stress  evolution  equation 

The  evolution  equation  for  Rij  is  derived  by  multiplying  the  ith  fluctuation  momentum  equation  by  u' 
and  the  j*1*  equation  by  u',  adding  the  resulting  equations,  and  averaging.  The  result  may  be  written  as 


fiRij  ~  Pij  +  0%)  d-  Tij  —  D{j  —  Jijktk  • 


(2.7.1) 


Here  the  production  term 

Pij  =  -(RikSkj  +  RjkSki)  (2.7.2) 

is  the  source  of  Reynolds  stress;  note  that  its  trace  is 

Pa  =  2  P.  (2.7.3) 

The  kinematic  rotation  term 

On  =  RikOkj  +  RjkOki  (2.7.4) 

is  trace  free  (O,,-  =  0)  and  hence  this  term  does  not  contribute  to  production  of  new  turbulence  energy,  but 
simply  rotates  the  turbulence  structure.  The  pressure  strain  term 


Ta  =  +-p'(u!.y  +“>» ) 

is  also  trace-free  and  provides  intercomponent  energy  transfer.  The  dissipation  term 


(2.7.5) 


Dij  =  2i/u',t  u'.,fc 


(2.7.6) 


has  a  trace 


Dij  =  IP. 


(2.7.7) 
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Finally,  the  flux  of  R,j  in  the  k,h  direction  i« 


Jijk  =  -(p'u'/.k  +  P'u'ijk)  +  u[u'u;  -  l/Rij,k  ■ 


(2.7.8) 


Again  we  have  used  the  mean  strain-rate  and  rotation-rate  instead  of  just  the  mean  velocity  gradients 
to  bring  out  the  different  roles  played  by  strain  and  rotation.  Most  previous  workers  have  included  the  mean 
rotation  term  in  with  the  production.  But  it  is  trace-free  it  does  not  add  new  energy  (it  is  absent  from  the 
turbulent  kinetic  energy  equation),  and  therefore  is  different  than  production.  The  rotation  term  provides 
exactly  the  changes  that  would  be  seen  if  the  Rij  structure  were  to  be  rotated  as  a  solid  body  without 
change.  Only  strain,  acting  on  the  Reynolds  stresses,  can  act  as  a  new  source  for  additional  Reynolds  stress. 

The  Rij  equation  is  sometimes  rewritten  so  that  the  trace  of  the  dissipation  term  is  2e  instead  of  20. 
with  an  associated  modification  in  the  flux.  This  is  left  as  an  exercise. 

The  R,j  evolution  equation  forms  the  basis  for  many  of  the  current  types  of  turbulence  models.  It  is 
very  useful  in  exploring  the  general  nature  of  the  changes  that  occur  in  turbulent  flows  subjected  to  strain. 


2.8  Vorticity  equation 

The  mean-square  turbulent  vorticity,  sometimes  called  the  enstrophy,  is  sn  important  quantity  in  tur¬ 
bulence.  Its  evolution  equation,  derived  by  multiplying  (2.5.2);  by  oi{  and  averaging,  is  useful  in  studying 
the  small  scales  of  turbulence.  Denoting 

(2.8.1) 


one  finds 


— ,1 


d(-u2)  =  uyjSij  -  u2Sjj  -  iv'u'n,,,  +n,s'yu>; 


•Hv.'w's'j  -  + 


(2.8.2) 


3.  STATISTICAL  DESCRIPTIONS  OF  HOMOGENEOUS  TURBULENCE 
3.1  Introduction 

A  field  in  which  all  statistical  properties  are  independent  of  position  is  called  homogeneous.  If  the 
statistical  properties  are  independent  of  the  orientation  of  the  coordinate  frame,  the  field  is  called  isotropic. 
Turbulence  may  be  approximated  as  homogeneous  and/or  isotropic,  although  turbulence  is  usually  homo¬ 
geneous  if  it  is  isotropic.  Few  practical  flows  are  either  homogeneous  or  isotropic.  Nevertheless,  regions  of 
practical  flows  often  are  essentially  homogeneous,  and  homogeneous  flows  provide  a  very  important  point  of 
departure  for  models  and  theories  of  turbulence.  Therefore,  development  of  good  understanding  of  homoge¬ 
neous  turbulence  is  an  important  first  step  in  the  study  of  turbulence. 

In  order  for  the  turbulence  to  be  homogeneous,  the  terms  in  the  equations  for  statistical  properties 
of  turbulence  must  be  independent  of  position,  Since  the  production  term  (2.7.2)  involves  mean  velocity 
gradients,  these  must  be  independent  of  position  if  homogeneity  is  to  be  achieved.  Therefore,  a  necessary 
condition  for  homogeneity  is  that  the  mean  velocity  be  a  linear  function  of  the  coordinates.  Since  there  are 
no  Reynolds  stress  gradients  in  homogeneous  turbulence,  the  mean  momentum  equation  (2.3.6)  shows  that 
the  mean  velocity  field  is  unaffected  by  the  turbulence. 

Since  the  mean  field  is  decoupled  from  the  turbulence  in  homogeneous  turbulence,  almost  any  mean 
velocity  history  can  be  imposed  on  homogeneous  turbulence.  Any  mean  Btrain-rate  history  can  be  imposed, 
but  the  mean  rotation  history  is  determined  by  the  imposed  strain-rate  history.  From  (1.10.5)  it  follows  that 
the  last  term  on  the  right  in  the  mean  vorticity  equation  (2.4.2)  is  (w'u(),y  which  vanishes  in  homogeneous 
turbulence.  Hence,  the  mean  vorticity  equation  in  homogeneous  turbulence  is 

n.  =nJs,J-n,s„,  (3.1.1) 

Thus,  while  an  initial  arbitrary  mean  rotation  can  be  imposed,  any  subsequent  changes  in  the  mean  rotation 
are  governed  by  (3.1.1).  This  restriction  is  important  in  the  analysis  and  simulation  of  turbulence  distortion 
by  mean  strain  and  rotation. 

The  statistics  of  homogeneous  turbulence  will  depend  upon  time.  Experiments  on  homogeneous  tur¬ 
bulence  generally  involve  passing  flow  through  a  grid,  which  generates  turbulence,  and  then  through  a  flow 
passage  of  varying  cross  section.  The  flow  is  approximately  homogeneous  as  seen  by  an  observer  moving 
downstream  with  the  mean  flow,  and  the  evolution  of  this  turbulence  as  seen  by  the  observer  is  interpreted 
as  the  time  evolution  of  the  turbulence.  The  behavior  of  turbulence  under  imposed  mean  strain  can  be 
studied  by  changing  the  cross-sectional  geometry  of  the  flow  channel.  Ingeneous  experiments  permit  great 
flexibility  in  such  experiments  (Gence  and  Mathieu  1980).  Homogeneous  shear  flow,  in  which  the  mean 
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streamwise  velocity  gradient  is  uniform  across  the  flow,  can  be  produced  using  upstream  grids  of  special 
geometry  (Tavoularis  and  Corrsin  1981). 

Grid  turbulence  is  not  quite  isotropic.  However,  by  placing  a  contraction  in  the  flow  stream  downstream 
of  the  grid,  essentially  isotropic  turbulence  can  be  obtained  (Comte-Bellot  and  Corrsin  1966)  for  study  in 
a  subsequent  constant  cross-section  duct  downstream.  One  can  also  study  the  relaxation  of  homogeneous 
turbulence  after  strain  in  such  a  duct. 

In  the  period  1960-1980,  a  number  of  basic  experiments  on  homogeneous  turbulence  provided  a  sound 
data  base  on  these  flows.  Since  then  numerical  simulations  of  homogeneous  turbulence  have  added  consid¬ 
erably  to  this  data  base.  The  insight  gleaned  from  these  experiment  and  simulations  now  allow  us  to  paint 
a  useful  picture  of  the  structure  of  homogeneous  turbulence.  The  next  section  presents  this  picture  and 
discusses  the  important  scales  of  turbulent  motion. 

3.2  Structure  and  scales  in  homogeneous  turbulence 

One  can  think  of  homogeneous  turbulence  as  a  complex  tangle  of  vortex  filaments,  each  acting  as  a 
“Biot-Savart  source”  in  moving,  distorting,  and  and  straining  all  the  filaments  (F:g.3.2.1).  This  continual 
vortex  stretching  concentrates  the  vorticity,  and  so  the  volume  of  vortical  fluid  tends  to  be  a  small  fraction  of 
the  total.  Vortex  filaments  of  the  same  sign  tend  to  collect,  and  this  provides  a  mechanism  for  the  creation 
of  larger  eddies.  This  is  counterbalanced  by  the  three-dimensional  straining  of  filaments,  which  tend  to 
twist  and  tangle  themselves  to  produce  smaller  eddies.  The  imposition  of  mean  strain  distort  the  tangle  of 
vortex  filaments,  much  as  the  fibres  in  a  ball  of  steel-wool  are  distored  when  it  is  stretched.  This  alters  the 
structure  of  the  filaments,  and  hence  the  structure  of  the  turbulence.  Upon  removal  of  the  mean  strain-rate, 
interactions  between  filaments  randomize  their  orientation,  bring  about  a  return  to  isotropy. 

This  tangle  of  vorticity  produces  a  very  broad  range  of  turbulent  motions.  The  larger  scales  are  more 
efficient  in  generating  the  Reynolds  stress  required  to  extranet  energy  from  the  mean  field  flow,  and  and 
new  turbulent  energy  appears  initially  at  large  scales.  Through  the  complex  non-linear  interactions,  which 
are  inviscid  processes,  turbulence  energy  is  cascaded  successively  to  smaller  and  smaller  eddies,  ultimately 
to  be  dissipated  by  viscous  straining  in  the  smallest  eddies,  where  the  local  strain  rates  are  the  greatest. 


Figure  3.2.1  Homogeneous  turbulence  as  a  tangle  of  vortex  filaments 


The  scale  of  the  largest  eddies  is  set  by  whatever  object  produced  them.  In  grid  turbulence  the  grid 
mesh  determines  the  largest  eddies,  in  wakes  the  large  eddies  scale  on  the  diameter  of  the  object,  and  in 
pipes  they  scale  on  the  pipe  diameter.  The  scale  of  the  smallest  eddies  is  set  by  the  rate  at  which  they  must 
dissipate  energy,  provided  to  them  by  the  large  eddies  through  the  cascade,  through  viscous  stresses.  The 
role  of  viscosity  in  turbulence  is  to  set  the  scale  of  the  smallest  eddies. 

These  ideas  suggest  that  the  dissipation  rate  is  determined  by  the  scale  of  the  energetic  large-scale 
turbulence  which  starts  the  energy  cascade.  If  we  assume  chat  q 3  and  e  characterize  these  large  scales,  then 
by  dimensional  analysis  the  length  scale  of  these  eddies  is 

l  =  q3/e  (3.2.1  a) 


and  the  time  scale  is 

r  =  ?3/e. 


(3.2.16) 


The  velocity  scale  is  of  course  just  q.  These  scales  tell  us  how  the  statistical  properties  of  large  eddies  should 
be  non-dimension alised  to  collapse  data  from  similar  flows  at  different  scales. 

The  Reynolds  number  of  the  turbulence,  defined  in  terms  of  the  velocity  and  length  scales  for  the  large 
eddies,  is 


(3.2.2) 


In  practical  flows,  q  is  generally  proportional  to  the  velocity  difference  driving  the  flow  (velocity  at  the  center 
of  a  pipe  or  the  velocity  defect  in  a  wake),  and  l  is  proportional  to  the  object  dimension.  Thus,  Rp  is  usually 
proportional  to  (but  smaller  than  by  a  factor  of  20-100)  the  flow  Reynolds  number. 
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The  scales  of  the  smallest  eddies  are  determined  by  the  e  and  v.  By  dimensional  analysis,  the  length 
scale  must  be 

Ik  =  (t'3/*)‘/4  (3.2.3a) 

and  the  time  scale 

tk  =  (^/e)1/3.  (3.2.3  6) 

These  Kolmogorov  scales  characterize  the  vortex  filaments  in  turbulence,  with  the  cores  of  the  vortex  being 
of  order  Ik  and  rotation  time  for  the  core  scaling  on  tk.  The  corresponding  velocity  scale,  characterizing 
the  velocity  difference  developed  locally  around  a  vortex  filament,  is 

VK  =  M‘/4 

Using  these  scales,  the  ratio  of  the  largest  scales  to  the  smallest  scales  is 

~  =  4/4  (32. .4) 

Thus,  the  range  of  turbulence  eddies  broadens  as  the  Reynolds  number  increases.  This  wide  range  limits 
direct  numerical  simulations  of  turbulence  to  low  Reynolds  numbers.  Large  eddy  simulations  of  turbulence, 
in  which  turbulence  of  smaller  scale  than  the  computational  mesh  is  modeled  and  the  larger  scales  are 
computed,  depends  heavily  on  models  for  the  small  scales.  It  is  tempting  to  approximate  this  sub-grid  scale 
turbulence  as  homogeneous,  and  therefore  a  firm  understanding  of  homogeneous  turbulence  is  important  to 
progress  in  large  eddy  simulation. 

The  remainder  of  this  chapter  is  devoted  to  the  mathematics  used  to  describe  the  statistical  properties  of 
homogeneous  turbulence.  Subsequent  chapters  deal  with  the  dynamic  evolution  of  these  statistical  properties 
in  response  to  imposed  mean  strain. 

3.3  Correlations  and  spectra 

The  statistical  properties  of  homogeneous  random  fields  are  most  often  described  in  terms  of  correlations 
and  spectra,  for  example,  if  /  and  g  are  two  random  field  variables,  the  two-point  correlation  of  /  and  g  is 

defined  as  _ 

Q/c (*.*'. 0  =<  f[x,t)g(x',t)  >  (3.3.1a) 

where  the  overline  denotes  a  volume  average  and  the  brackets  denote  an  ensemble  average.  Ensemble  and 
volume  averages  are  usually  assumed  to  be  the  same  for  homogeneous  fields  ( ergodic  hypothesis);  the  dual 
averaging  is  therefore  redundant  but  useful  in  the  analysis  that  follows. 

For  homogeneous  fields  Qf,,  depends  only  on  the  separation  of  the  two  points  r  =•  (x'  -  x)  and  t, 

Q/y{r >0  =<  f(x,t)g(x  +  r,t)  >  .  (3.3.16) 

Often  the  time  dependence  of  the  correlation  is  not  expressed  explicitly,  but  it  must  not  be  forgotten. 

There  is  an  infinite  set  of  other  correlations  of  possible  interest,  for  example  the  two-point  correlation 
with  time  delay,  three-point  correlations,  etc.  A  complete  statistical  description  requires  knowledge  of  the 
probability  density  function  for  all  variables  of  interest  at  all  space  points  and  time,  ar.  impossible  goal  to 
achieve.  Therefore,  statistical  descriptions  are  always  limited  in  what  they  can  provide,  and  the  challenge  is 
to  provide  what  is  really  essential,  with  minimum  effort  and  maximum  accuracy. 

In  homogeneous  fields,  Fourier  expansions  can  be  used  to  represent  individual  realizations  of  the  fields. 
Suppose  that  /  and  g  are  defined  within  a  box  of  interest  (Fig.  3.3.1).  In  order  to  give  them  Fouriei 
expansions  we  have  to  imagine  that  they  are  periodic  functions,  so  let 

/(x)  inside  the  box 

periodic  repetition  outside. 

The  Fourier  representation  of  /  at  any  instant  of  time  is 

/(x)  =  £/>')*-“■' x 

k' 


(3.3  2a) 


(3.2.3c) 


where  k  =  (tj,^,^)  is  the  three-dimensional  wavenumber  vector,  and  k  x  =  knxn.  Since  the  Fourier 
modes  must  fit  into  the  box  with  integer  periods, 

ki  =  2irnj/L.  (3.3.26) 

The  summation  is  a  triple  sum  over  all  Fourier  modes, 
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£  =  £££ 

k'  k  i  kj  kt 

Note  that  the  Fourier  coefficient*  may  vary  with  time;  we  do  not  ehow  this  explicitly  here. 


(3.3.2c) 


Figure  3.3.1  Box  for  Fourier  representation 

There  is  an  important  relationship  between  the  Fourier  coefficients  of  positive  and  negative  wavenumbers 
for  a  real  field.  Taking  the  complex  conjugate  of  (3.3.2a),  replacing  k'  by  k", 

/•(*)  =  £/*(  k'Vil,Hx 

k" 

where  the  *  denotes  a  complex  conjugate.  Letting  k"  =  -k', 

/*(*)  =  £/*(- k,)e-V*.  (3.3.3) 

k' 

Now,  if  /  is  real  it  is  equal  to  its  complex  conjugate.  Equating  the  Fourier  coefficients  of  like  exponentials 
in  (3.3.2a)  and(3.3.2), 

/>')  =  /*(- k') 

or  al‘  ernatively  (for  real  f ) 

/(-k)  =  /*(k).  (3.3.4) 

The  Fourier  coefficients  are  evaluated  using  the  orthogonality  property  of  the  Fourier  modes,  using 
integrals  over  the  domain.  In  what  follows  f()d3x  denotes  an  integral  over  the  box  in  Fig.  3.3.1.  Then, 
multiplying  (3.3.2a)  by  e’kx  and  integrating  over  the  box, 


J  e,k  ’7(x)d3x  =  £/(k')  J  e,<k-k'>  Xrf3x. 


Since  each  Fourier  mode  that  fits  the  box  contains  an  integer  number  of  cycles, 


/  .(k-k')  xjs  fO  ifk^k' 

J  \  L3  ifk  =  k'. 


Hence  all  terms  in  the  summation  of  (3.3.4)  drop  out  except  for  the  one  where  k'  =  k.  Thus,  the  Fourier 
coefficients  can  be  evaluated  as 

/( k)  =  jj]  f(x)eikxd3x.  (3.3.7) 

The  two-point  correlation  of  /  and  g  can  be  expressed  in  terms  of  the  Fourier  representations.  Consider 
the  correlation  of  /  and  g  within  the  box  of  Fig.  (3.3.1), 

<  /(x)j(x')  >=  ££  <  /(k)j(k')  >  e-«k  x+k'  x'). 


The  brackets  indicate  that  the  Fourier  coefficients  are  random  variables  that  will  differ  from  from  realisation 
to  realisation.  Let  k"  -  -k'  and  r  =  x'  -  x.  Then,  using  (3.3.4), 

<  /(x)j(x  +  r)  >=  ££  <  /(k)?*(k")  >  e-<x  <k-k">+>k"  (3.3.8) 


Now  we  average  by  integrating  over  the  box  and  dividing  by  L3,  denoting  this  average  by  an  overline.  Using 
(3.3.6),  all  the  terms  in  the  sum  drop  out  except  the  terms  where  k"  =  k.  The  result  is  then 
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$/.( r)  =<  /» j(x  +  r)  >=  £  <  />)5*(k)  >  «*k  r-  (3.3.9) 

k 

In  computational  simulations  in  which  the  evolution  of  the  Fourier  coefficients  is  calculated  for  finite-series 
approximations  to  the  fields,  (3.3.9)  is  used  to  calculate  the  the  two-point  correlation. 

Theoretical  treatments  take  the  limit  as  L  ->  oo,  in  which  case  the  sums  become  integrals.  To  pass  to 
this  limit,  we  note  that  the  difference  between  consecutive  wavenumbers  in  the  summation  is  Afca  =  2ir /L 
for  each  direction,  so  Afc,L/2ir  =  1.  We  can  multiply  each  term  in  the  summation  by  unity  three  times  to 
obtain 

Qfg[r)  -■=  E  <  /(k)$‘00  >  (£)  UciMiAkz'*  '.  (3.3.10) 

We  define  the  cospectrum  of  /  and  g  as 

®/',k,=  (^)S</"(k)r(lt)>'  (3311) 

This  is  the  equation  used  to  calculate  the  co-spectrum  in  discrete  spectral  simulations  of  random  fields. 
Then, 

0/,(r)  =  £  E/9(k)eik  'Afc,Ak2A*3.  (3.3.12) 

k 

Taking  the  limit  as  L  — *  oo,  we  define  the  cospectrum  of  /  and  g  by 

£/„(k)  =  Urn  £/(l(k).  (3.3.13) 


E/g  does  not  become  infinite  as  L  — •  oo  because  the  Fourier  coefficients  of  individual  modes  become  very 
small  as  the  number  of  significant  modes  increases.  As  L  — •  oo,  AfciAfc3Afc3  becomes  an  elemental  volume 
in  wavenumber  space  dkidkjdkj  =  d3k.  Therefore,  in  (3.3.12)  the  two-point  correlation 

Q/Ar)  =  ^  <?/<,(«•)  (3-314) 

becomes  the  three- dimensional  Fourier  transform  of  the  cospectrum, 

<?/.( r)  =  /  F/,(k)e<k  *<f3k.  (3.3.15) 

Here  the  triple  integration  is  to  be  carried  out  over  all  wavenumbers. 

There  is  an  inverse  of  the  transform  (3.3.15).  Multiplying  (3.3.9)  by  e~,k  r  and  integrating  over  a  box 
of  size  L  in  r  space, 

I  $„(r)«-v- rd3r  =  E/  <  /WW  >  «ir  ^k■k',  (3.3.16) 

Each  exponential  in  the  summation  will  execute  an  integer  number  of  cycles  in  each  direction  and  hence 
integrate  to  zero,  except  for  the  term  where  k  =  k'.  Hence, 

/  Q/»(r)«-<k'  r<*3r  =  L3<  /( k)j-(k)  >=  (2*)3£/*(k).  (3.3.17) 

In  computational  simulations  based  on  finite-difference  methods,  this  equation  is  used  to  calculate  the  cospec¬ 
trum  from  the  two-point  correlation.  Taking  the  limit  as  L  ->  oo,  and  replacing  k'  by  k, 


k)  = 


(3.3.18) 


Note  that.  E ' j g  and  Qjg  are  Fourier  tranform  pairs. 

We  could  have  obtained  the  cospectrum  simply  by  Fourier  transformation  of  the  two-point  correlation. 
We  started  with  a  finite  box  so  that  the  relationships  between  the  Fourier  coefficients  and  the  cospectrum 
would  be  made  clear,  and  also  to  derive  results  useful  to  persons  engaged  in  discrete-representation  simula¬ 
tions  of  homogeneous  turbulence  in  finite  computational  domains.  It  should  be  understood  that  the  Fourier 
transforms  of  /  and  g  defined  over  an  infinite  region  do  not  exist.  However,  because  events  at  distant 
separations  are  uncorrelated,  Qfe  — ►  0  as  |r|  — *  oo,  and  hence  the  Fourier  transform  of  Q/s  does  exist. 
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з. 4  Velocity  correlations  and  apeetra 

The  velocity  field  in  homogeneous  turbulence  can  be  represented  in  the  terms  outlined  above,  Let  /  = 

и,  and  g  =  uy.  Then,  dropping  the  redundant  ensemble  average, 


<?■>(«•)  =  u'(x)u;.(x  +  r).  (3.4.1) 

Qij  is  the  two-point  velocity  correlation  tensor.  Note  that 

9«(0)  =  u;(*K(x)  =  q2.  (3.4.2) 

Qij(f)  expresses  the  average  relationship  between  two  velocity  components  measured  at  two  locations 
separated  by  a  distance  x.  Qaa  ( repeated  Greek  indices  are  not  summed)  will  be  largest  for  sero  separation, 
fall  to  a  fraction  of  this  value  for  separations  comparable  with  the  large  eddies  in  the  turbulence,  and  become 
sero  for  infinite  separation.  If  the  eddies  tend  to  be  long  in  one  direction  and  short  in  another,  this  will 
be  reflected  in  the  different  rate  at  which  the  correlation  falls  off  with  different  ra.  Thus,  the  two-point 
correlation  tensor  can  tell  one  quite  a  bit  about  the  structure  of  the  turbulence. 

Using  (3.3.18),  the  velocity  spectrum  tensor  is 

£.,(k)=  (^)7  «.>)<-’- rd3r  (3.4.3) 

where  the  integrations  are  over  all  r.  It  is  related  to  the  two-point  velocity  correlation  tensor  by  (3.3.15), 

<?.y(r)  -  j  fty( k)e<k  rd»k  (3.4.4) 


where  the  integrations  are  over  ail  h. 

The  Reynolds  stresses  ft,  =  u'u'  are  given  by 

ft,  =  <?.,(o) 


(3.4.5) 


for  which  (3.4.4)  gives 

ft,  =  J  fty(k)d3k.  (3.4.6) 

Reviewing  the  developments  of  the  previous  section,  one  sees  that  ft,(k)d3k  represents  the  contributions  to 
fty  coming  from  an  element  of  k  space  of  volume  d3 k  positioned  at  k. 

For  uniform  density  flow,  the  continuity  equation  (2.3.4)  provides  important  constraints  on  Q,y(r).  from 
(3.4.1) 

=  u'(x)u'(x  + r),y  =  0.  (3.4.7a) 

Replacing  x  by  x*  —  r  in  (3.4.1),  then  differentiating  with  respect  to  r,-,  (2.3.4)  also  requires  that 


(3.4.76) 


The  continuity  equation  (2.3.4)  also  constrains  fty.  In  terms  of  the  Fourier  expansion,  continuity 
»-equires 

-5>yfly(k)«"<fc-“  =  o.  (3.4.8) 

k 

This  must  hold  for  all  x,  which  requires  that  the  coefficient  of  each  and  every  exponential  must  vanish. 
Hence,  for  each  wavenumber  vector  k, 

fcyiiy(k)  =  0.  (3.4.9) 

Equation  (3.4.9)  is  the  continuity  equation  in  Fourier  form.  It  says  that,  for  each  k,  the  Fourier  coefficient 
vector  <1  must  be  orthogonal  to  k  in  order  for  the  velocity  field  to  be  divergence-free.  This  condition  is  used 
very  often  in  analysis  and  simulation  of  homogeneous  turbulence.  From  (3.4.9),  it  follows  (most  obviously 
using  the  the  discrete  Fourier  representations)  that 


KEij  =  0 


(3.4  10a) 
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and 

kjEi,=  0.  (3.4, 106) 

The  correlation  tensor  Qi,  has  an  important  symmetry  property.  Noting  that 

Q.j(-r)  =  o' (x)u' (x  -  r) 

we  let  x  =  x'  +  r  and  rewrite  this  as 

Q>j(-r)  =  u'(x'  +  r)u'(x'). 

The  right  hand  side  is  just  Q,,  (r).  Hence 

=  QM*)  (3  4.H) 

The  spectrum  tensor  Eij  also  has  a  symmetry  property.  Since  the  Fourier  coefficients  for  real  fields 
obey  (3.3.3),  It  follows  (most  obviously  from  the  discrete  Fourier  representations)  that 

4(-k)  =  (^:)9<fli(-kK{-k)> 

=  (^)3  <  fii(k)Mk)  >“  (3.4.12) 

In  the  limit  L  — *  oo  this  becomes 

£«(-k)  =  ^(k).  (3.4.13) 

The  turbulence  kinetic  energy  may  be  expressed  as 

\<?  =  \q..(0)  =  \  I  Bui k)d3k.  (3.4.14) 

Integral  scales  of  motion  may  be  defined  in  terms  of  Q*j-  For  example, 


rQn(ri,0,0)dr, 

“  <?n(0,0,0) 


(3.4.15) 


is  useful  as  a  measure  of  the  Xi  scale  of  the  turbulence.  Here  the  arguments  display  the  three  components 
of  the  separation  vector, 

Qij  and  Eij  are  the  classical  quantities  used  to  describe  homogeneous  turbulence.  They  are  less  use¬ 
ful  for  inhomogeneous  turbulence  because  expansion  functions  other  than  Fourier  modes  really  should  be 
used  in  directions  of  inhomogeneity.  They  are  used  for  inhomogeneous  Rows  when  the  turbulence  can  be 
approximated  as  locally  homogeneous  over  regions  large  compared  to  the  integral  scale. 


3.S  Other  statistical  quantities 

There  are  many  other  statistical  quantities  of  interest  in  turbulence.  Those  that  involve  only  quadratic 
forma  in  the  velocity  are  termed  second-order.  Any  second-order  statistical  property  of  the  velocity  field  can 
be  derived  from  the  two-point  correlation  tensor  or  the  velocity  spectrum  tensor.  For  example,  a  tensor  of 
interest  is 

^•JPS  =  uJtpuyn-  (3.5.1) 


PYom  (3.4.1), 


dQ.Ar) 

dr. 


(x  +  r). 


(3.5.2) 


Replacing  x  by  x'  -  r  in  (i  5.2),  then  differentiating  with  respect  to  rp,  one  has 


drpr, 


-u',p(x'-r)U;„  (x'). 


(3.5.3) 


drpdr,  )  |r|_0 


(3.5.4) 


Now  letting  r  =  0, 


1-19 


The  corresponding  result  in  terms  of  the  spectrum  tensor  can  be  derived  directly  by  taking  the  derivatives 
of  the  discrete  Fourier  series  for  the  velocities,  and  proceding  as  in  section  3.2  above,  or  by  applying  (3.5.4) 
to  (3.4.4).  Either  approach  gives 

Dijpq  =  J  kpk,Ei, (k)d3k.  (3.5.5) 

Since  gradients  of  all  statistical  quantities  vanish  in  homogeneous  turbulence, 

K“y)»y=°-  (3.5.6) 

Expanding  the  differentiation  using  the  continuity  equation  (2.3.4), 

Dim  =  “!■>“>..  =  0-  (3.5  7) 

Note  that  this  is  consistent  with  (3.5.4)  and  (3.5.5)  if  the  continuity  constraints  (3.4.7)  or  (3.4.10)  are  applied. 
The  dissipation  rate  e  may  be  expressed  in  general  as 

e  =  u(Diiii  +  Diiii).  (3.5.8) 

FYom  (3.5.7),  the  second  term  does  not  contribute,  and  in  homogeneous  turbulence  the  true  dissipation  rate 
t  is  the  same  as  the  homogeneous  dissipation  rate  P  defined  by  (2.6.4). 

Using  (3.5.8),  (3.5.7),  and  (3.5.5),  we  find  that  the  dissipation  rate  is  related  to  the  velocity  spectrum 
tensor  by 

t  =  u  j  *2£„  (k)d3k.  (3.5.9) 

The  factor  k2  means  that  the  main  contributions  to  the  dissipation  come  from  higher  wavenumbers  (smaller 
eddies)  than  those  that  provide  the  major  contribution  to  the  kinetic  energy. 

8.6  Vorticity 

The  two-point  vorticity  correlation  tensor  is 


Note  that 


FYom  the  definition  of  vorticity  (1.9.1), 


H'.jfr)  =w'(x)u»'(x  +  r) 


W..(0)  =  =  wa. 


U  —  Dim  —  Dt 


so  it  follows  from  (3.5.7)  and  (3.5.8)  that  in  homogeneous  turbulence  the  dissipation  is  directly  related  to 
the  mean-square  vorticity, 

e  =  vw2 .  (3.6.4) 

The  enstrophy  equation  (2.8.1)  is  therefore  sometimes  used  as  a  guide  in  developing  model  equations  for  the 
dissipation. 

The  vorticity  can  also  be  expanded  in  a  Fourier  representation;  for  the  box  of  section  3.2, 


=  £  w** 


Because  the  vorticity  is  by  definition  divergence-free. 


and  because  the  vorticity  is  real 


k,w,  (k)  =  0 


<*00  =  w*(-k). 


The  vortxcxty  spectrum  tensor  /f,y(k)  can  be  developed  following  the  approach  above.  It  is  of  course 
the  Fourier  transform  of  the  two-point  vorticity  correlation  tensor,  and  can  be  related  to  the  velocity  tensor. 
Because  the  vorticity  is  divergence-free, 

kiHijl  k)  =  0  (3.6.8a) 


*/ff«(k)  =  0 


(3.6.85) 


1-20 


and  because  it  is  real 

JM— k)  =  %(k). 

The  Fourier  coefficients  of  the  vorticity  are  related  to  those  of  the  velocity.  Using  (1.9.1), 

*'(x)  =  x. 


(3.6.10) 


Equating  coefficients  of  like  exponentials,  the  vorticity  coefficients  are  found  to  be 

ut'(k)  =  ~*kp£,-pgtg(k). 


(3.6.11) 


From  this  it  follows  that 


f f ,j(  k  )  —  r  i  kpkr  Eq ,  (k  ) 


(36.12) 


One  can  express  <?,y  in  terms  of  the  vorticity  correlation  tensor  and  E,y  in  terms  of  the  vorticity  spectrum 
tensor.  This  requires  the  solution  of  the  Poisson  equation  (1.9.5),  which  is  easily  accomplished  using  the 
Fourier  representations.  Alternatively,  one  can  multiply  (3.6.10)  by  er,ik,.  The  result  is 


“i(k)  =  fip,-rf(i,(k). 


(3.6.13) 


where  ft2  =  Jyky.  Substituting  in  the  discrete  representation  of  EtJ  and  taking  the  limit,  one  finds 

E,y(k)  =  f,.p,£yr.^tf„(k). 


(3.6.14) 


This  result  finds  important  use  in  rapid  distortion  theory,  where  it  is  used  to  estimate  the  anisotropy  in  the 
Reynolds  stresses  produced  by  distortion  of  the  vorticity  field  due  to  imposed  mean  strain.  It  is  also  useful 
in  constructing  models  of  E, y  for  anisotropic  turbulence  using  models  for  the  anisotropic  //,y. 

3.7  Correlations  and  spectra  in  isotropic  turbulence 

If  the  statistics  are  independent  of  the  coordinate  system  orientation,  only  two  types  of  correlations 
completely  characterise  the  velocity  correlation  tensor  (Fig.  3.7.1).  The  longitudinal  correlation  function 

/(r)  =  4<3*‘(r*-°-°)  (3.7.’) 

r 

describes  the  coherence  of  the  velocity  fluctuations  aligned  with  the  separation  of  the  two  points.  The  lateral 
correlation  function 

g(r)  =  -~Q22(n>0,0)  (3.7.2) 

<T 

relates  to  the  coherence  of  fluctuation  velocities  perpendicular  to  the  separation  axis. 


Figure  3.7.1  Longitudinal  and  lateral  correlation  functions 

The  complete  tensor  Qiy(r)  can  expressed  in  terms  of  these  two  scalar  functions.  The  tensor  must  be  a 
function  of  the  separation  vector  r.  The  most  general  such  function  is 

Q.y(r)  =  Cx6ij  +  C3r,ry  (3.7.3) 

where  the  coefficients  C\  and  Cj  may  be  functions  of  the  scalar  invariant  of  the  vector,  r  =  s/rfr\.  The 
coefficients  can  be  identified  by  expressing  the  longitudinal  and  lateral  correlations: 


Q»(ri,0,0)  =  \g(r)  =  C1 


(3.7.4) 
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Qi  l  (rl 1 0, 0)  =  j  /(r)  =  Ci  +  Cjr3, 


Solving  for  Ci  and  Cj,  one  finds 


QiAr)  = 


/(')  ~  g(r) 


ur,  +  g(r)fij 


(3.7.5) 


(3  7.6) 


Note  that  /  and  g  are  scalar  functions  of  the  scalar  separation  magnitude  r  (and  of  time,  not  shown  explicitly). 
The  continuity  equation  provides  a  relationship  between  /  and  g.  Since  r  =  i/firi, 


dr  r* 
dr*  r  ' 


Differentiating  (3.7.6)  with  repect  to  r*, 


Qii*  =  J  i~^)  r,r->  7'  +  (~r^)  +  +  9‘~ 


(3.7.6) 


(3.7.8) 


where  the  primes  denotes  differentiation  with  respect  to  r.  Setting  k  =  j  and  using  the  continuity  condi- 
tion(3.4  7),  one  finds 

/'+!(/- 9)  =  0  (3.7.9) 


This  integrates  readily  to  give 


fir)  =  4  f  rg(f)df. 
r  Jo 


(3.7.10) 


Eij  for  isotropic  turbulence  can  be  obtained  by  Fourier  transform  of  Qt)  as  outlined  in  section  3.3. 
Alternatively,  we  can  construct  its  general  form  directly  since,  for  isotropic  turbulence  the  Eij  tensor  must 
be  a  function  only  the  vector  k.  The  most  general  form  is 


£„(k)  =  CiS.j  +  Cjk,k,  (3.7.11) 

where  the  coefficients  can  depend  on  the  scalar  invariant  of  the  vector,  k  =  \/(kiki)  Using  the  continuity 
condition  (3.4.10), 

Cikj  +  C7k7kj  =  0  (3.7.12) 

hence 

Co  =  -Ci/k2  (3.7.13) 

Redefining  C\  as  4ir k2E{k),  we  have 

<»•»> 

E(k)  is  called  the  energy  spectrum  [unction.  Note  that  it  is  a  scalar  function  of  the  scalar  k  (and  of  time, 
not  shown  explicitly). 


Figure  3.7.2  Coordinate  system  for  fc-space  integration 
The  turbulence  energy  is,  using  (3.4.14), 

b*  =  /  <3-715) 


The  integration  of  integrals  of  this  type  ,  in  which  the  unknown  function  depends  only  on  the  magnitude  of 
the  vector,  can  best  be  carried  out  in  spherical  coordinates  (Fig.  3.7.2).  We  have 
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i^±r  r  r  mk,tin$dedHk. 

*  4 n  Jk=0  Jo  =  Q  J 

Carrying  out  the  integrations  over  <f>  and  9, 


=  j~  E(k)dk. 


(3.7.16) 


(3.7.17) 


We  see  that  E(k)dk  represents  the  contribution  to  the  kinetic  energy  per  unit  mass  arising  from  all  the 
Fourier  inodes  in  a  spherical  shell  in  k-space  of  radius  k  and  thickness  dk.  Once  E(k)  is  known,  the  entire 
velocity  spectrum  tensor  Etl  is  known  from  (3.7.14). 

In  theory,  homogeneous  isotropic  turbulence  evolves  in  time,  and  one  should  measure  the  spectrum 
tensor  by  making  measurements  at  many  space  points.  In  reality  this  is  very  difficult  (but  it  is  what  is 
in  fact  done  in  a  numerical  simulation).  Instead,  laboratory  experiments  make  use  of  Taylor’s  hypothesis, 
which  assumes  that  the  velocity  pattern  measured  as  a  function  of  time  at  one  point  is  frozen  in  the  fluid 
and  being  swept  over  the  probe.  The  probe  measurement  is  thereby  interpreted  as  providing  Qn(ri,0,0). 
Using  (3.7.14)  in  (3.4.4), 

Qn(riiO, 0)  =  /  ( 1  -  p)eifc‘r‘i3k  (3.7.18) 

This  integration  is  conveniently  carried  out  in  the  coordinates  of  Fig.  3.7.3.  We  sort  the  Fourier  contributions 
according  to  those  with  the  same  wavenumber  |Atf  | .  Terms  from  both  sides  of  the  kt  axis  contribute,  with 
opposite  signs  in  their  exponentials;  these  are  combined  into  a  cosine: 


(3.7.19) 


f3n(rli°i0)  =  fitm0fMi  lm0  S'  (l  -  |)2cos [klrl)kd*dkdkl 

We  carrying  out  the  <j>  integration,  and  define  the  one  dimensional  spectrum  function  Ei  by 

E>{ki)=L-h  (3-7-2o) 

Then, 

<3u(ri,0,0)  =  f  Ei[k)cos(kiri)dki  (3.7.21) 

Jk,=  o 

One  can  taking  the  Fourier  cosine  transformation  of  the  measured  <?i i (*"x ,  0, 0)  to  get  £"!  (A: j ) .  Then, 
differentiating  (3  7.18)  twice  (a  courageous  step  with  laboratory  data!), 


r.,,  ,  kla'Esjki)  k\  dE\  (<C|) 

[  l]~  2  dk\  2  dkt 


(3.7.22) 


This  allows  E(k)  to  be  determined.  It  also  shows  that  if  E(k)  varies  as  a  power  of  k  in  some  range  the.i 
Ei(ki)  will  vary  with  the  same  power  of  k\. 


Figure  3.7.3  Coordinates  for  one-dimensional  spectrum  integration 

Eqn.  (3.7.21)  in  essence  defines  Qu  as  a  one-dimensional  Fourier  cosine  transform  of  El.  The  inverse 
transform  is 

o  roo 

(3.7.23) 


Ei(ki)  =  -  [  <?n(ri,0,0)cos(*:iri)dri. 

*  Jo 


Noting  that  Qn(0)  =  q*/ 3  in  isotropic  turbulence,  the  integral  scale  defined  by  (3.4.15)  is  given  by  (3.7.23) 


in  -  A /  =  —Ei(0) 


(3.7.24) 
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A/  is  the  integral  scale  derived  from  the  longitudinal  correlation  function  /(r),  and  hence  it  is  called  the 
called  the  longitudinal  integral  scale.  Since  it  is  non-sero,  E\  (0)  >  0,  in  contrast  to  E( 0)  =  0. 

3.8  Diaaipation  Id  Isotropic  turbulence 

Using  the  isotropic  spectrum  (3.7.14)  in  (3.5.9),  carrying  out  the  integrations  using  polar  coordinates 
as  above,  the  dissipation  rate  is  found  as 


t  =  i/  f°°  k?E(k)dk 

Jo 


The  factor  k2  means  that  higher  wave  numbers  (smaller  scales)  make  more  contribution  to  the  dissipation 
(and  vorticity)  than  they  do  to  the  energy  (compare  3.7.17). 

Since  the  small-scale  component  of  turbulence  is  generally  throught  to  be  very  nearly  isotropic  at  high 
Reynolds  numbers,  isotropic  turbulence  theory  is  used  as  an  aid  in  estimating  e  from  laboratory  data.  This 
approach  makes  use  of  the  tensor  Z),/p«  defined  by  (3.5.1).  In  an  isotropic  field,  the  only  tensors  upon  which 
DiJpq  can  depend  are  the  isotropic  numerical  tensors,  hence  it  must  be  of  the  form 

Dijpq  —  CiSjjSpq  +  CjSipSjq  +  C3S,qSIP  (3.8.2) 

where  the  coefficients  must  be  scalars.  The  coefficients  can  be  evaluated  from  three  known  constraints.  First, 

the  definition  forces  a  symmetry, 

Ajm  =  Awp*  (3.8.3a) 

Second,  continuity  requires  that 

Dijiq  =  0.  (3.8.36) 


Finally,  we  know  that 


Using  these  considition,  one  finds 


e  =  vD« 


(3.8.3c) 


2  £  1 

Aim  =  j  [^oA>«  ~  J . 


The  pertinent  quantity  most  easily  measured  in  an  experiment  (again  using  Taylor’s  hypothesis)  is 

GCT-lfc’  <3-8'5) 

This  is  usually  the  way  that  e  is  estimated  in  laboratory  experiments. 

Another  important  turblence  scale  defined  in  terms  of  the  dissipation  is  the  microscale.  It  can  be 
approached  through  the  longitudinal  correlation  function  /(r).  The  symmetry  property  (3.4.11)  indicates 
that  /(r)  must  be  an  even  function  of  r,  so  its  expansion  is 

/(r)  =  1  -  ~ar2  +  0(r4)  (3.8.6) 

The  interception  of  this  osculating  parabola  (Fig.  3.8.1)  with  /  =  0  defines  a  scale  A;  =  y/2/a,  called  the 
longitudinal  Taylor  microscale.  From  (3.5.4),  using  (3.7.5)  and  then  (3.8.4), 

3  n  e 

0  -  -jt-'im  -  ~ — n 
q*  5i jq* 

so 

X)  =  lOi rq2/e.  (3.8.7) 

Alternatively,  the  dissipation  rate  can  be  expressed  as 

q 2 

e  =  10i/-rj.  (3.8.8) 

This  equation  is  sometimes  used  to  determine  e  from  measurements  of  the  longitudinal  correlation. 

Using  (3.2.3)  and  (3.2.2),  the  ratio  of  the  Taylor  microscale  to  the  Kolmogorov  scale  is 


^  =  yio4/4 

Ik 


(3.8.9) 
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Using  (3.2.1),  the  ratio  of  the  energy-containing  scale  to  the  Taylor  scale  is 

-L  =  -4=1#*  (3.8.10) 

A/  vTo  ^ 

Thus  the  microscale  falls  between  the  smallest  and  largest  scales.  Although  it  is  the  most  commonly  reported 
turbulence  scale,  it  is  the  least  well  understood.  It  has  been  suggested  that  it  is  a  measure  of  the  size  of  the 
loops  in  the  vortex  filaments,  but  this  is  not  at  all  certain. 


Figure  3.8.1  The  osculating  parabola  defines  the  Taylor  microscale 

8.9  Scaling  of  the  spectrum  in  isotropic  turbulence 

The  general  form  of  E(k)  deduced  from  measurements  in  isotropic  turbulence  is  shown  in  Fig.  3.7.1.  By 
(3.7.17),  the  area  under  the  curve  is  the  turbulent  kinetic  energy,  to  which  then  greatest  contributions  come 
from  wavenumbers  around  the  peak.  The  vorticity  and  dissipation  occur  predominantly  at  high  wavenumbers. 


m 


k 


Figure  3.9.1  Form  of  the  spectrum  in  isotropic  turbulence 

It  is  generally  thought  that  the  small-scale  motions  in  any  turbulent  flow  become  isotropic  at  high 
Reynolds  numbers,  and  therefore  that  the  Kolmogorov  scales  characterize  the  high  wavenumber  region 
of  any  turbulent  flow.  Moreover,  if  one  assumes  that  there  is  a  universal  small-scale  spectrum,  then  by 
dimensional  analysis  it  must  be  of  the  form 


The  one-dimensional  spectrum  £i(fci)  would  have  to  scale  in  the  same  manner.  Figure  3.9.2  shows  that  the 
data  from  a  wide  variety  of  flows  do  indeed  collapse  when  plotted  in  these  Kolmogorov  variables.  The  data 
flatten  at  low  wavenumbers  because  they  are  one-dimensional  spectra  where  £i(0)  is  given  by  (3.7.24). 
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STREAMWISE  ENERGY  SPECTRA  FOR  VARIOUS  TURBULENT  FLOWS 

(CHAPMAN.  1979) 
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Figure  3.9.2  Spectra  in  Kolmogorov  variables 


Kolmogorov  suggested  that  there  should  be  a  range  of  wavenumbers  in  which  the  main  process  is  the 
passing  of  energy  from  larger  eddies  to  smaller  eddies  (the  cascade  of  turbulence  energy),  and  that  the 
structure  of  this  region  should  depend  only  on  the  rate  of  energy  cascade.  Since  this  cascade  ultimately  ends 
with  dissipation  at  the  small  scales,  the  rate  of  energy  cascade  must  be  e.  If  one  assumes  that  E(k)  depends 


only  on  t  (and  of  count  k)  in  this  range,  by  dimenaional  analyaia 


J5(fc)fcs^3e  7^3  =  constant  =  a 


or 

E(Jfc)  =  ae3/3k~&/3.  (3.9.2) 

This  is  the  Komogorov  spectrum,  a  cornerstone  of  turbulence.  Measurements  give  a  Kolmogorov  constant  a 
of  about  1.5.  The  data  of  Figure  3.9.2  show  the  -5/3  range,  with  longer  runs  of  -5/3  behavior  at  larger 
Reynolds  numbers,  consistent  with  the  broadening  of  scales  as  . 

In  the  vicinity  of  the  peak  in  E(k),  the  spectrum  should  scale  on  the  large-scale  variables  (see  section 
3.1),  and  hence  should  collapse  when  plotted  as 


(3.9.3) 


Where  this  form  overlaps  with  the  Kolmogorov  spectrum  the  function  G  must  be  such  that  q  drops  out,  and 
this  again  establishes  the  -5/3  law  for  the  asymptotic  overlap  range  between  low  and  high  wavenumbers. 

Figure  3.9.1  indicates  that  E(k)  ->  0  as  k  -*  0,  but  there  is  controversy  as  to  just  how.  There  are 
good  arguments  supporting  both  k3  ( Saffman )  and  k*  ( Loitskianski )  variation  as  k  0.  The  kA  behavior 
is  required  if  £?,•/  is  to  be  analytic  at  k  =  0.  The  k3  behavior  implies  some  residual  preferential  directions 
at  lero  wavenumbers,  which  may  be  more  characteristic  of  physical  experiments.  Numerical  simulations 
with  delta  spectra  at  mid-range  611-out  as  it 4  as  the  turbulence  develops,  but  simulations  initiated  with  k2 
behavior  at  low  wavenumbers  persist  as  k3.  Simple  models  of  turbulence  show  that  the  rate  of  energy  decay 
in  isotropic  turbulence  depends  on  the  low  wavenumber  portion  of  the  spectrum,  and  with  the  experimental 
decay  rates  support  the  k3 . 


k 

Figure  3.9.3  Evolution  of  the  spectrum  in  decaying  isotropic  turbulence 

Turbulence  not  subjected  to  mean  deformation  will  decay  as  time  passes.  The  larger  eddies  take  more 
time  to  change,  and  the  smallest  scales  of  motion  adjust  most  rapidly.  Figure  3.9.3  depicts  the  nature  of 
the  evolution  of  E(k,t)  (we  now  include  the  time  variable  heretofore  suppressed).  Note  that  the  peak  moves 
to  larger  scales  (small  wavenumbers)  because  the  smaller  eddies  die  out  faster.  Thus  as  time  progresses  the 
integral  scale  will  grow. 


Figure  3.9.4  Model  spectrum  for  isotropic  turbulence 

A  simple  model  spectrum  for  isotropic  turbulence  is  shown  in  Fig.  3.9.4.  It  assumes  a  power  law 
behavior  at  low  wavenumbers,  a  -5/3  inertial  range,  and  a  sharp  cutoff  at  the  Kolmogorov  scale: 

(  Akm  for  k  <  ki, 

E(k)  =  |  ac2l3k~bl3  for  kL  <  k  <  kv  .  (3.9.4) 

(  0  for  k  >  kv 


Matching  the  spectrum  at  kL  gives 


-(=?) 


i/3\  3/(3m  +  S) 


Assuming  fc„  >■  kL, 


\q7  =  L  mdk  =  a{-^Ti  +  i)*i2/v/’ 


from  which  an  estimate  of  the  peak  wavenumber  is  obtained, 


wr2  /  i  +n]3/a. 

e  \m+  1  2/ 


Again  assuming  kv  3>  fct,  the  viscous  cutoff  wavenumber  is  estimated 


from  which 


=  v  f  k2 E{k)dk  =  v-ac2l3k*l 
Jo  4 

Kv3!*  _  (  4  \3/* 

«*/<  \3a/ 


It  is  left  as  an  exercise  to  work  out  the  one-dimensional  spectrum  E\  for  this  model  spectrum,  and  from  that 
to  determine  the  integral  scale.  For  m  =  2  and  a  =  1.5,  one  finds 

A/®/?3  “  0.11.  (3.9.9) 

This  model  spectrum  exhibits  the  proper  scaling  for  isoropic  turbulence,  and  gives  values  of  the  scales  within 
about  a  factor  of  two  of  those  found  from  actual  spectra.  It  is  very  useful  in  constructing  simple  turbulence 
models,  in  setting  up  initial  fields  for  turbulence  simulations,  and  in  addressing  other  aspects  of  homogeneous 
turbulence. 
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3.10  Third-order  statistics  in  isotropic  turbulence 

Correlation*  involving  product*  of  three  quantitie*  are  called  third  order  statistics.  These  depend  on  the 
relative  phase*  of  the  Fourier  modes,  information  not  contained  in  the  spectrum  tensor  E%).  Of  particular 
interest  in  turbulence  modeling  are  one-point  third-order  statistics.  For  isotropic  turbulence  these  tensors 
can  be  worked  out  using  the  methods  used  previously.  For  example, 

Zyfi,  =  CqWtiik.  (3.10.1) 

In  dealing  with  the  vorticity  and  dissipation  equations,  one  encounters  the  tensor 

d’ijkpqr  =  uJip  uyn  u*ir-  (3.10.2) 

This  is  evaluated  for  isotropic  turbulence  by  first  writing  the  general  tensor 

V'i/fcpfr  =  lip{Cl&je&kr  +  C3Sjl,Sqr  +  C3S]rSkq)  +  Sij(CtSpq  6  hr  +  C36pySqr  +  CaSprSqt) 
+^e(Oj6PJ6i,r  +  Ct6pi,6jr  +  C3Spr6jk)  +  ^,k(C,0SPJ  Sqr  +  Cn6pq6lr  +  Cl36pr6iq) 


+&\r(Cis6pj6qk  +  Cu6pq6jk  +  CiiSpkSqj)' 

(3.10.3) 

There  are  three  are  three  symmetry  constraints, 

fyjkpqr  ~  V* jipqkr 

(3.10.4a) 

V* ijkptjr  =  $kjirqp 

(3.10.46) 

V* ijkpqr  =  $ikjprq’ 

(3.10.4c) 

Continuity  also  provides  some  constraints,  but  with  the  symmetries  enforced  on'v  one  is  required, 

VUj  kipq  —  0. 

(3.10.5) 

Forming  («'«’ and  using  homogeneity  conditions,  one  can  show  that 

=  0. 

(3.10.6) 

With  these  constraints,  (3.10.3)  can  be  reduced  to  a  form  containing  only  one 
Cj  =  A,  one  finds 

unknown  coefficient.  With 

^Pijkpqr  —  ^  [{^ip^jq^  kr  *4*  ^ij^pk^qr  *4"  ^ij^pr^qk  ~4  ^iq^pr^jk  ^ik^pj^qr  "4“  ^ik^pq^jr  "4*  ^ir^pq^jk) 

4  3  1 

T  &qr  “f"  &i  j &pq &k r  "4"  &ik  ^pr  &jq )  —  T  ( ^iq ^pk  &jr  “4"  &ir$p]&qk  )  —  T  &kq  "4"  &iq 

a  4  0 

Spjhkr  +d.rdpk6j.i)]-  (3.10.7) 

For  example, 

K..  )3  =  * 

(3.10.8) 

u'u'  t'  ■  =  -  —  A. 

t  j  ij  2 

(3.10.9) 

The  derivative  tkewneu  is  related  to  A;  using  (3.8.5)  and  (3.10.8), 

tK,.]’/K7F1,/1-4(^)  '. 

(3.10.10) 

The  skewness  is  measured  to  be  negative,  the  term  given  by  (3.10.9)  is  positive.  This  is  the  turbulent  vortex 
stretching  source  term  in  the  equation  for  mean-square  vorticity  (2.8.2),  by  which  the  turbulence  tends  to 
enhance  its  own  mean-square  vorticity. 


4.  RAPID  DISTORTION  OF  HOMOGENEOUS  TURBULENCE 


4.1  Introduction 

The  state  of  homogeneous  turbulence  changes  significantly  when  it  is  subjected  to  mean  strain.  This 
occurs  in  practice  whenever  turbulence  passes  through  a  duct  of  variable  cross-section,  such  as  a  nozzle,  when 
turbulence  is  sheared  by  the  mean  flow,  or  when  turbulence  is  subjected  to  a  mean  rotation,  The  general 
trends  can  be  understood  using  vortex  stretching  concepts.  For  example,  passing  turbulence  through  an 
axisymmetric  nozzle  stretches  the  vortex  filaments  in  the  flow  direction  and  tends  to  align  them  with  the 
flow  direction,  reducing  the  turbulent  fluctuations  in  the  direction  of  flow  but  increasing  the  fluctuations 
transverse  to  the  flow. 

Because  of  the  non-linearity  of  the  governing  equations,  it  is  impossible  to  develop  a  rigorous  theory  of 
these  processes.  There  are  two  alternative  approaches  to  such  a  theory.  The  first  is  to  use  some  sort  of  a 
closure  model  to  produce  a  set  of  closed  equations  describing  the  evolution  of  statistical  properties  of  the 
turbulence  in  response  to  the  mean  strain.  The  second  approach  is  to  simplify  and  then  solve  the  exact 
equations  for  special  cases.  Both  approaches  are  useful.  In  this  chapter  we  examine  rapid  distortion  theory 
(RDT),  in  which  the  exact  equations  for  the  fluctuation  field  are  approximated  in  a  way  that  is  valid  for 
very  strong  imposed  mean  strain  rate,  yielding  linear  equations  amenable  to  exact  solution. 

It  might  be  thought  that  the  response  to  large  strain  rate  could  be  calculated  using  the  Reynolds  stress 
transport  equations  (2.7.1),  neglecting  the  terms  that  do  not  explicitly  contain  the  mean  velocity  gradients. 
However,  this  analysis  overpredicts  the  changes  in  the  Reynolds  stresses,  because  the  pressure-strain  term 
Tij  in  (2.7.1)  produces  an  immediate  effect  that  reduces  the  impact  of  the  production  term  P by  a  factor 
of  about  50%.  The  Poisson  equation  for  the  fluctuation  pressure  (2.5.3)  shows  that  a  sudden  onset  of  mean 
velocity  gradient  instantly  changes  the  fluctuation  pressure  field.  The  result  is  a  sudden  change  of  Ttj  with 
the  onset  of  applied  Siy,  and  this  must  be  considered  in  the  analysis.  Turbulence  modelers  refer  to  the  part 
of  T,j  that  changes  suddenly  with  a  sudden  change  in  the  mean  deformation  rate  as  the  rapid  pressure  strain 
term.  RDT  plays  a  key  role  in  understanding  and  modeling  this  term,  and  this  chapter  is  intended  to  aid 
the  use  of  RDT  in  this  work. 

The  basic  idea  of  RDT  is  that  when  |5|q2/e  >  1  the  time  scale  of  the  turbulence  q2 /e  is  long  compared 
to  that  of  the  mean  deformation,  and  so  the  turbulence  does  not  have  time  to  interact  with  itself.  Thus,  the 
non-linear  terms  in  the  governing  equations  (2.5 . l)-(2.5.3)  involving  products  of  fluctuation  quantities  are 
neglected,  and  so  the  RDT  equations  are  linear  in  the  fluctuation  quantities.  The  viscous  terms  are  linear 
and  can  be  included  in  the  analysis,  but  are  often  neglected  and  will  be  here. 

These  equations  contain  the  mean  velocity  gradients,  which  must  be  independent  of  position  for  ho¬ 
mogeneous  turbulence  but  may  depend  on  time.  The  convective  operators  D  contain  the  mean  velocities, 
which  must  vary  linearly  with  x  in  homogeneous  turbulence.  These  coefficients  prevent  representation  of 
the  solution  as  periodic  in  the  coordinates,  and  this  hampers  direct  solution  by  Fourier  methods.  However, 
when  transformed  to  coordinates  marked  on  the  mean  flow  at  the  start  of  the  distortion,  the  variable  co¬ 
efficients  are  removed  and  the  solution  may  be  obtained  by  Fourier  methods  in  the  transformed  system. 
This  transformation  is  used  in  the  numerical  simulations  of  homogeneous  turbulence  (Rogallo  1981),  where 
it  permits  the  numerical  solution  to  be  ezacf  for  infinitely  rapid  distortions!  The  numerical  simulations  of 
the  full  equations  carried  out  using  this  program  are  useful  in  helping  assess  the  range  of  applicability  of 
RDT,  and  it  is  rather  surprising  that,  for  some  types  of  strain,  RDT  works  remarkably  well  even  at  relatively 
low  strain-rates  (Lee  and  Reynolds  1985).  Thus,  RDT  is  becoming  recognized  as  being  very  important  and 
useful  in  turbulence  analysis,  modeling  and  simulation  (Savill  1987). 

4.2  The  RDT  equations 

The  most  general  mean  velocity  field  in  which  homogeneous  turbulence  can  exist  is  of  the  form 

U,  =  Aik(t)xk  (4.2.1a) 

from  which 

Ui,k=A,k(t).  (4.2.16) 

Note  that  (3.1.1)  restricts  the  rotational  history  of  the  imposed  mean  deformation,  but  any  mean  strain  rate 
history  can  be  imposed. 


I  'll 


Substituting  (4.2.1)  in  (2.5. t ) ,  the  inviscid  RDT  approximation  for  homogeneous  turbulence  is 

u'  +  Ajkxku'itj  =  -u'  A,j  -  -p'„  .  (4  2.2) 

P 

The  continuity  equation  (2.3.4)  also  applies.  We  remind  the  reader  that  these  equations  hold  if  p  =  p(t),  so 
they  can  be  applied  in  certain  types  of  compressible  flow  situations. 

Solution  by  Fourier  methods  is  practical  only  if  the  coefficients  in  the  equations  are  independent  of  x. 
Therefore,  it  is  necessary  to  transform  the  equations  to  remove  the  troublesome  term  on  the  left-hand  side. 
The  transformation  is  assumed  to  be 


(i  =  Bik{t)ik  T-t. 


(4.2.3) 


TVansforming  (4.2.2)  to  the  new  coordinates,  the  left-hand  side  becomes 

1 9u,  dui  ■  dui  „ 

17  +  irnB"kXk  +  A,kXkduBnj' 

Setting  the  coefficient  of  x*  to  zero  to  remove  the  variable  coefficient, 


Bnk  +  AjkBnj  =  0. 


(4.2.4) 


This  defines  the  Rogallo  transformation.  The  B,j  can  be  found  by  solving  these  linear  equations,  although 
a  closed-form  solution  is  not  feasible.  The  transformation  ties  the  new  coordinate  systems  to  the  mean 
motion,  with  the  new  grid  distorting  and  rotating  as  demanded  by  the  mean  flow.  The  Rogallo  code  for 
direct  simulation  of  homogeneous  turbulence  operates  in  this  coordinate  system. 

With  this  transformation  the  RDT  momentum  equations  (4.2.2)  become 


£4 

dr 


-u'd„ 


dp 


and  the  continuity  equation  (2.3.4)  becomes 


3u[ 


d(k 


-Bki  =  0. 


(4.2.5a) 


(4.2.56) 


The  Poisson  equation  for  the  pressure  fluctuation  is  obtained  by  taking  the  derivative  of  (4.2.5a)  with  respect 
to  and  the  derivative  of  (4.2.5b)  with  respect  to  r  and  combining,  using  (4.2.4).  Alternatively,  one  can 
simply  transform  (2.5.3).  The  result  is 


1  3V  3u' 

--zrir8"8™  =  (4.2.5c) 

These  linear  equations  can  be  solved  to  track  the  evolution  of  the  Fourier  coefficients  of  the  velocity  field  in  the 
transformed  coordiates.  The  Reynolds  stresses  are  integrals  of  this  spectrum  function,  and  the  integrations 
may  be  cariied  out  in  the  transformed  coordinates.  If  the  spectrum  in  the  original  coordinates  is  involved, 
the  spectrum  must  be  mapped  back  to  x  space  using  the  coordinate  transformation. 

Closed-form  solution  of  the  RDT  equations  for  a  general  problem  is  not  possible.  However,  exact 
solutions  for  special  cases  can  be  obtained,  in  some  cases  in  closed  form  and  in  others  in  terms  of  integrals. 
The  general  solution  can  be  found  as  a  power  series  in  time.  Some  of  these  solutions  that  play  useful  roles 
in  understanding  turbulence  and  in  turbulence  modeling  will  now  be  discussed. 

4.3  Response  of  turbulence  to  rapid  rotation 

RDT  can  be  applied  to  study  the  effect  of  rapid  rotation  on  turbulence  in  the  absence  of  strain.  Taking 
the  rotation  as  clockwise  about  the  13  axis,  the  mean  velocity  is 


Ui  =  Txj 


Uj  =  -r*i 


(4.3.16) 


and  the  coordinate  transformation  is  (Fig,  4.3.1) 


£ i  =  ii  cos(rr)  -  ijsinfTr) 

(4.3.2a) 

(3  =  ijcos(rr)  +  11  sin(rr) 

(4.3.26) 

(3  =  13 

(4.3.2c) 

r  =  (. 

(4.3. 2d) 

Transforming,  the  RDT  equations  become 

17  =  ^[cos(rT)^  +  8in(rr)^]“U3r 

(4.3.3a) 

^  in(rr)^-  +  cos(rr)^]  +  Ulr 

(4.3.36) 

duj  1  dp 

dr  pd£3 

(4.3.3c) 

The  transformed  continuity  equation  is 


du\  /„  ,  3u'.  .  ._  .  duU  .  ,  3u',  du'3 

cos(rr)  +  —  sm(rr)  -  —  ..»(!>)  +  —  cos(rr)  +  ^  =  0. 


(4.3.3d) 


It  is  helpful  to  transform  the  velocity  components  to  the  rotating  coordinate  system.  Denoting  these 
velocities  by 

Vi  =  ui  cosfTr)  -  u2sin(rr)  (4.3.4a) 

v2  =  u2cos(rr)  +  uj  sin(Tr)  (4.3.46) 

u3  =  u3.  (4.3.4c) 

Forming  the  equations  for  the  new  velocities  form  the  old,  one  finds 


T1 r 

dr  pdt  i 

dv2  1  dp  „  _ 

7T  =  —pW2+1 

du3  1  dp 
dr  pd(3 

dv>  „ 


(4.3.5a) 

(4.3.56) 
(4.3.5c) 
(4.3. 5d) 


The  second  terms  on  the  right  are  of  course  the  Coriolis  terms. 

Now  we  seek  the  solution  for  the  evolution  of  the  Fourier  modes  in  the  transformed  space.  Following 
the  developments  of  section  3.3,  we  write 


K 


(4.3.6a) 


P  =  (4.3.66) 

K 

where  k  is  the  wavenumber  in  the  transformed  coordinates  Equating  coefficients  of  like  exponentials, 

diii  \K\  . 

—  =  — p-  2Tv2 
dr  p 


(4.3.7a) 
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dOj  ‘*3  .  .  ,r. 

__  =  — p  +  2rv! 
or  p 

(4.3.76) 

303  t  /c3  „ 

_  =  — p 

or  p 

(4.3.7c) 

* 

e> 

II 

p 

(4.3.7d) 

Applying  the  continuity  equation  (4.3.7d)  to  (4.3.7a-c), 

1  .  2r(02/c,  -  fiiicj) 

P  4 

p  t*4 

(4.3.8) 

where  /c3  =  /c3  +  x3  +  x3.  Substituting  (4.3.8)  in  (4.3.7a-c),  and  seeking  solutions  of  the  form  fh(x,  r)  = 
a,exp(i/Jr),  one  obtains 

—  2F  “(a2Xi  —  fljx2)  +  2Ta2  =  0  (4.3.9a) 

\pa.i  -  2r^|(a2X[  -  a|X2)  -  2Tai  =  0 

(4.3.96) 

t/9a3  -  2r^|(a2/Ci  -  a,K2)  =  0. 

X* 

(4.3.9c) 

This  linear  equation  system  has  non-trivial  solutions  only  if  the  determinant  of  the  coefficient  matrix  vanishes. 
This  condition  gives 

/?3  =  4r3^f  -  j  =  4r3^|  >  0.  (4.3.10) 

Hence,  except  for  modes  with  x3  =  0,  the  solutions  are  undamped  oscillations  in  time  at  frequency  /?(x). 

The  x3  =  0  modes  require  special  attention.  They  correspond  to  two-dimensional  modes  with  their 
vorticity  aligned  with  the  rotation  axis.  The  solution  for  these  modes  is 

<'i(s.r)  =  SiU.O)  ~  C(s)K2r 

(4.3.11a) 

Ms,  0  =  Ms,  °)  +  <?(ts)«i  r 

(4.3.116) 

>>3(5,  f)  =  v3(«,  0) 

(4.3.11c) 

where 

c  „r^«lCi('S,0) +K3fi2('S,0)  j 

(4.3. lid) 

But  for  x3  =  0  the  numerator  of  C  is  lero  by  continuity,  and  hence  the  Fourier  coefficients  of  these  modes 
do  not  change  under  rapid  rotation.  Thus,  these  coefficients  can  also  be  regarded  as  undamped  oscillations 
at  frequency  /J(x). 

The  solution  for  the  Fourier  coefficients  is  therefore 

v,  =  a+e’'^+a.e-’^. 

(4.3.12) 

aj±  and  a2±  are  related  by  (4.3.9a)  or  (4.3.9b), 

\±l~  +  — )«1±  =  (^3  “ 

(4.3.13) 

The  coefficients  a,±  are  set  by  the  initial  values  of  the  Fourier  amplitudes, 

t>i0  =  u«+  +  ®v- 

(4.3.14) 

where  0,o  is  the  initial  value  of  G,(x).  Using  (4.3.13)  and  (4.3.14),  one  finds 

Ql±  =  [(T’7  +  ^)c,° +  (!  ■ 

(4.3  15) 

Following  section  3.4,  the  spectrum  tensor  £;,  (in  the  rotating  coordinate  system)  is 


£.,(*. r)  =  <  ft  (a.  »■)*/(«.»■)  >  • 


(4.3. 1C) 


Using  the  solution  and  a  bit  of  algebra,  one  finds 


£u('S,r)  =  +  ^4~ )  £n(«,0)  +  2^1  -  j  £22(1,0)  +  ^1  -  ^^(£12(1,0)  +  £21(1,0)) 

+2[(^  "  -  (l  -  ^|)^22(l,0)  -  ^(1  -  §)(£.2(l,0)  +  £21(1,0))]  cos(2/?r) 

-2[^i^£,i(ic,0)  +  ^(1-^)(£i2(l,0)  +  £2i(l,0))  sin(2/?r)  j.  (4.3.17) 


-2^^p£,i(/c,0)  +  -f  [1  -  (£12(1,0)  +  £21  (l,0)) J  ain(2/?r)  j.  (4.3.17) 

If  the  initial  turbulence  is  isotropic,  the  initial  spectrum  is  given  by  (3.7.14),  and  one  finds  that  the 
coefficients  of  the  sin  and  cos  terms  vanish;  hence  there  is  no  change  in  the  spectrum  as  viewed  by  an  obsever 
in  the  rotating  coordinate  system.  Since  the  spectrum  is  isotropic,  the  spectrum  seen  by  a  stationary  observer 
is  also  unchanged.  Thus,  rotation  of  itself  will  not  distort  the  spectrum  of  isotropic  turbulence. 

If  the  initial  spectrum  is  anisotropic,  as  for  example  produced  by  prior  strain  and  associated  rotation, 
the  residual  rotation  will  simply  cause  the  spectrum  to  oscillate  at  a  frequency  w  —  2/3(/c).  The  associated 
Reynolds  stress  (in  the  rotating  frame),  determined  by  integrating  £n  over  all  5,  will  oscillate  in  a  compli¬ 
cated  manner  that  depends  on  the  initial  spectrum.  However,  using  the  symmetry  property  of  the  spectrum 
(3.4.13),  the  contribution  of  the  sin(2/3r)  term  to  the  integral  is  seen  to  vanish.  Hence,  relative  to  a  rotating 
observer,  the  Reynolds  stress  oscillations  can  be  expressed  as  an  even  power  series  in  r  arising  from  the 
cos(2 fir)  term.  Hence,  the  Reynolds  stresses  seen  by  a  stationary  observer  would,  to  O(t),  appear  to  rotate 
in  the  manner  described  by  the  kinematic  rotation  terms,  with  deviations  from  this  behavior  being  described 
by  an  even  power  series  in  time. 

These  are  important  results  for  turbulence  modeling.  Turbulence  models,  when  reduced  to  the  same 
rapid  distortion  approximations,  should  not  show  any  effect  of  pure  rotation  (rotation  without  straining) 
on  isotropic  turbulence.  Moreover,  when  applied  to  the  pure  rotation  of  anisotropic  turbulence,  the  models 
should  shown  the  kinematic  rotation  of  the  Reynolds  stress  described  by  (2.7.4),  plus  modifications  by  an 
even  power  series  in  time.  This  condition  is  very  useful  in  setting  coefficients  in  turbulence  models,  and  we 
shall  use  it  in  Chapter  6. 

4.4  Rapid  isotropic  compression  or  expansion 

Consider  next  isotropic  expansion  (or  compression)  with 


The  RDT  momentum  equations  are 


U,  =  fz,. 


+  r xu'itj  =  -Tu'  -  ^yp'„  • 


The  density  is  given  by  the  continuity  equation, 


p  =  -3pl\ 


The  RDT  transformation  is 


and  the  transformed  equations  are 


£,  -  z,«' 


d-A  =  -lV- 
dr  1  p(i)a£,' 
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*4 

<?£, 


=  o. 


(4.4.7) 


Multiplying  (4.4.6)  by  u|,  the  pressure  term  dropa  out  by  continuity  (4.4.7).  Averaging,  we  obtain  the 
RDT  approximation  for  the  kinetic  energy, 


2  dr  H 


(4.4.8) 


The  aolution  ia 

=  qlr 7rt  (4.4.9) 

where  q 1  ia  the  initial  kinetic  energy.  Thua,  the  turbulence  kinetic  energy  will  decrease  with  expansion 
(T  >  0)  and  increase  with  compression. 

The  evolution  of  the  spectrum  ia  obtained  by  solving  the  individual  component  equations.  Fourier 
expansions  are  used  as  above.  The  pressure  fluctuations  (i.e.  the  rapid  part)  are  zero  by  continuity,  and  all 
Fourier  modes  of  the  velocity  vary  as  ex p(-Tt).  Thus,  the  spectrum  retains  its  inital  shape  in  the  etretcked 
coordinate  system,  and  simply  scales  in  magnitude  with  q 2.  As  a  consequence  the  integral  scale  (3.4.15) 
varies  in  proportion  to  the  strain, 

A/(0  =  A/(0)er‘.  (4.4.10) 

These  results  are  useful  in  constructing  turbulence  models  for  compressible  turbulence.  Some  of  the 
turbulence  models  currently  in  use  do  not  predict  the  proper  behavior  with  compression,  some  even  predicting 
an  increase  in  length  scale  as  turbulence  is  compressed! 

4.5  Response  of  turbulence  to  rapid  Irrotational  strain 

RDT  analysis  for  irrotational  mean  strain  is  neatly  handled  using  the  vorticity  equation.  Under  the 
RDT  approximations,  with  no  mean  rotation,  (2.5.2)  reduces  to 


Dw[  =  w'S.j  -  w-Stn. 

We  work  in  principal  coordinates  of  Si)  and  take 

u*  =  r„(t)*.. 


(4.5.1) 


(4.5.2) 


Recall  that  Greek  indices  are  not  summed.  The  RDT  coordinate  transformation  is 


(a  —  Xa/ea  T  —  t 


(4.5.3a, 4) 


where 


is  the  total  strain  in  the  a  direction.  The  transformed  vorticity  equation  is 


(4.5.3c) 


where 


The  solution  of  (4.5.4)  is 


du'°  -  f 

TT  ~  rnU,« 

(4.5.4a) 

fa  =  To  —  To 

(4.5.5a) 

r0  =  Tr  +  r2  +  r3. 

(4.5.56) 

ua(X’T)  =  “'Llx.Ojea 

(4.5.6) 

a  =  exp  fa(t')dt'^ 

(4.5.7) 

where 


I-.VS 


is  a  modified  total  strain  in  the  a  direction.  This  result  clearly  shows  the  essence  of  RDT;  it  computes 
the  change  in  the  turbulent  state  by  considering  the  rapid  vortex  atreching  imposed  by  the  mean  field.  The 
velocity  field  can  be  deduced  from  the  vorticity  field.  In  the  transformed  coordinates,  the  Poisson  equation 
(1.9.5)  for  the  velocity  gives 

- -j  ,  dui,  1  1  . ,  .  _  . 


W,  =  ^ 

1  Q  c  -  AC 


where  the  transformed  Laplace  operator  is 


V2  = 


«3  d^2  e2 


Id2  Id2  Id2 


t  ^  „  .o  "f  n  ~ 


t\  t\ail  e2  d{2 


(4  5.8a) 


(4.5.86) 


The  equations  for  U)'2  and  w'3  can  be  obtained  by  permuting  the  indices.  The  solution  is  obtained  using 
Fourier  expansions, 

w'( x,  r)  =  ^c5j(s.,  r)e_,*"("  (4.5.9a) 


«!(*, T)  =  “•(*.  r)e""c”(’‘ 


(4.5.96) 


The  solution  for  fii  is 


j(/c3w2/e3  -  x2w3/e2) 

2  2  2‘ 

72  +  +  -y 


(4.5.10) 


The  other  components  can  be  lound  by  permutation  of  the  indices. 

The  velocity  spectrum  function  F,y  is  related  to  the  vorticity  spectrum  function  Hi y  by 

r,  <  >  (^)5^22(x,0)  +  (i%)2 H33(h,0)  -  2(^)(^)H23(k,0) 

£,n(x,r) -  - -  . 


From  the  solution  for  the  vorticity  evolution  (4.6,6), 


Hu/»(x,  r)  =  HaP(K,Q)eaC0. 


(4.5.12) 


If  we  assume  that  the  initial  turbulence  is  isotropic,  the  initial  vorticity  spectrum  is  given  by  (3.7.23), 
with  k  replaced  by  x.  Using  this  spectrum  and  (4.5.12)  in  (4.5.11), 


£it(s.r) 


£(x)  (f?)M(*Z  ~  *1)  +  +  2(|ft)x!x2 


(4.5.13) 


The  spectra  of  £22  and  £33  can  be  found  by  permuting  the  indices. 

The  Reynolds  stresses  can  now  be  calculated  by  integrating  £*y  over  all  wavenumbers  (see  3.4.6).  The 
integrations  are  most  easily  carried  out  using  spherical  coordinates,  and  can  be  evaluated  in  closed  form  for 
a  few  very  simple  cases,  such  as  isotropic  compression.  However,  the  general  case  of  irrotational  strain  can 
be  handled  by  power  series  expansion  in  the  total  strains.  In  (4.5.3c)  we  expand 


<„  =  exp(a)  =  1  +  aa  +  -a2  +  . 


(4.5.14) 


The  integrals  are  then  expressed  as  power  series  in  the  aa,  and  evaluated  in  spherical  coordinates,  where  the 
angular  integrations  can  be  carried  out  analytically.  The  x  integration  produces  ?2/3,  the  initial  isotropic 
value  of  Rn.  Using  this  approach,  the  Reynolds  stress  Rjj,  dissipation  tensor  Z^y,  and  vorticity  Vn  =  wjcvl 


I  -  v> 


were  evaluated  by  the  author  and  Moon  Lee  (Reynolds  1983)  to  0(a2).  Subsequently  Piomelli  used  the 
symbolic  manipulation  program  MACSYMA  to  extend  the  Reynolds  stress  (the  most  important  quantity 
for  turbulence  modeling)  to  0(a5)  (Lee,  Piomelli,  and  Reynolds  1986).  The  results  are  as  follows  (the  other 
components  can  be  found  by  permuting  the  indices): 


^11  —  ?o 


5"  i(2a°  +  4a‘ 


+  ttt  ( 10flo  +  12a2  -  32a2a3)  + 
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-  12a2  +  48aoa2Q3  —  8a 2  +  48a;a2a3) 
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225225' 


d~  -(30a*  +  32agOi  —  48a2a2  —  128a2a2a3  —  320aoa;a2a3  -+  36aJ  —  32a2a2a3  +  16a2a3) 

3465 

-244a';  -  1440ajai  +  3360ajX  -  4320ajjaJ  -  2488a5  +  1200a^a2a3  +  7600a(2)a1a2a3  +  4320aoa* 


-3440aoa2a2a3  -  880aoa2a2  +  5200a5a2a3  +  560aia|a2)  +  0(a°) 
fin  —  2co 


12  1  2, 

---aI  +  -a0  +  O(a) 


V,i  =  y 


e  =  eo 


1  +  2aj  —  ay  4-  O(a^) 


1“  jao  +  0(a2) 


»n  =  +  2aM)  +  ~K3  -  wJ)  +  -^K4  +  4a;2a*a;  -  2a;2a;2) 
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where  the  Reynolds  stress  anisotropy  tensor  is 


(196819a;5  +  25189a;a;2aJ2  -  479453a;3a;a;)  +  0(a'°) 


i  fit;  72^t;/3 

(y,  = - 5 — — 

72 


(4.5.15) 

(4.5.16) 

(4.5.17) 

(4.5.18) 

(4.5.19) 

(4.5.20) 


and  the  anisotropic  strain  components  are 


;  a»  -  ~"n- 


(4.5.21) 


Note  that  the  anisotropy  tensor  6,y  is  depends  only  on  the  total  anisotropic  strain,  and  is  independent  of  the 
strain-rate  history.  These  results  are  useful  in  turbulence  modeling  where  one  seeks  to  develop  models  that 
will  be  consistent  with  RUT  when  the  RDT  approximations  are  applied  to  the  model. 

4.6  Combinations  of  strain  and  rotation 

The  general  RDT  problem  for  homogeneous  turbulence  involves  combinations  of  strain  and  rotation,  for 
which  a  general  solution  can  be  developed  in  symbolic  form  (Cambon  1981).  Using  the  Fourier  expansions 
(4.5.9)  and  a  similar  one  for  the  pressure,  (4.2.5c)  is  first  solved  to  express  the  Fourier  coefficients  of  the 
pressure  in  terms  of  those  of  the  velocity, 


1.  2iKkBkiAij  . 

~P  =  ~  n - 5*-“; 

P  *-n*rr\  £}my>Dnp 

Then,  the  Fourier  expansion  of  (4.2.5a)  gives 


dii 

dr 


=  -AijUj  + 


2  KkBkiK.,,B,w 
Kn^m  Bm,Bn, 


ApjUj  —  HijUj. 


(4.6.1) 


(4.6.2) 


Following  Cambon,  the  solution  can  be  expressed  using  Green’s  functions, 


u;(«,r)  =  Ci  k(K,r)uk(K,0) 


(4.6.3) 


where  the  Green’s  functions  are  given  by  the  solution  of 

30  U,r)  //|t(gtf)c  (g[T)  o  'a, 0)  =  (4.6.4) 

OT 

This  allows  the  spectrum  tensor  Eij  to  be  expressed  in  terms  of  the  initial  spectrum, 

EuU,  r)  =  G\,,(/c,  r)G‘q(K,  t)Em(k,0).  (4.6.5) 

The  Reynolds  stresses  are  then  simply  integrals  of  the  spectrum  function  over  all  <j. 

This  method  of  solution  is  instructive  for  looking  at  the  structure  of  the  solution,  but  the  calculations  for 
the  Reynolds  stresses  require  approximate  evaluation  of  the  integrals,  for  example  by  power  series  expansions. 
Moreover,  when  the  principal  axes  of  the  strain  rate  vary  with  time  the  Green’s  functions  are  not  easily 
obtained,  except  perhaps  as  power  series  in  time.  If  one  is  going  to  resort  to  series  solution,  a  direct  solution 
by  power  series  in  time  is  simpler.  We  will  develop  this  here  for  future  reference. 

A  superscript  summation  convention  aids  the  analysis.  We  denote  a  series  by 

OO 

A  =  £  A[n)t"  =  A(n)tn.  (4.6.6) 

n=0 

Any  repeated  superscript  or  power  is  summed  over  all  possible  values.  The  delimiters  on  the  superscripts 
establish  the  range,  with  (  )  establishing  a  lowest  value  of  0  and  [  )  establishing  a  lowest  value  of  unity. 

Muliplication  of  two  series  and  sorting  out  of  powers  of  t  is  then  very  easily  accomplished.  For  example, 
simply  replacing  n  by  r  -  m  in  the  product  below  collects  the  coefficient  of  tr, 

AB  =  =  A[r'm)B{m)tr.  (4.6.7) 

Here  the  delimiters  correctly  establish  that  the  m  summation  in  the  coefficient  of  tr  is  from  0  to  r.  The 
leading  coefficients  can  also  be  extracted, 

A^-m)Btm)  =  ^(rlgfO)  +  At<>)  gtr)  +  ^Ir-m)  fl|.n|  (4.6.8) 

where  the  m  sum  at  the  end  now  ranges  from  1  to  r  —  1. 

We  treat  a  general  case  of  arbitrary  strain  and  initial  rotation  as  applied  to  initially  isotropic  turbulence, 
and  express  the  velocity  gradient  tensor  A,k  (see  4.2.1)  as 

<4.*(t)  =  £*(<)  +  ^£«,n;,(  t).  (4.6.9) 

The  strain-rate  history  described  by  5,*(t)  and  the  initial  rotation  described  by  0,  (0)  will  be  arbitrary,  and 
the  rotation  history  is  governed  by  (3.1.1). 

Expanding, 

sik  =  s^tn  n,  = 

Aik  =  /»!*'«"  Blk  =  sj^r 

U,'  =  p  =  p|n)r. 

The  coefficients  then  are  generated  recursively.  FYom  (3.1.1) 


(4.6.10) 


(4.6.12) 


l-.tx 


FVom  (4.6.11) 


FVom  (4.2.4) 


j('l  _  o(r)  .  If  fjM 


(r+i)o!;;l'  =  -/>'r,l<1  = «-*- 

The  Fourier  representation  (in  stretched  space)  of  (4.2.5a)  gives 
(r+  l)u,(r  +  11  = 


The  continuity  condition  (4.2.5b)  gives 

_.  .-(r-»)  n(v)  _  n 
K*ui  "t»  =  "■ 

The  Fourier  representation  of  the  Poisson  equation  (4.2.5c)  gives 

Extracting  the  leading  pressure  term  given  by  r  =  0, using  (4.6.15b), 

--xV  =  -*tKmP,r"')l  |rlBir*,Bm’  +  2t«Ct  *>  Ui>] 

P  P 


(4.6.13a,  fc) 

(4.6.14) 

(4.6.15) 

(4.6.16) 

(4.6.17) 


where  the  notation  llrl  forces  r  >  0  in  the  sum. 

The  procedure  is  now  very  simple.  At  each  order  r,  one  finds  the  rotation  term  from  from  (1.6.11),  the 
velocity  gradient  term  from  (4.6.12),  the  transformation  term  from  (4.6.13),  the  pressure  term  from  (4.6.17), 
and  the  velocity  terms  from  (4.6.14).  The  spectrum  tensor  is  expressed  as  a  similar  series  expansion,  and  its 
terms  are  generated  and  integrated  in  spherical  coordinates  to  calculate  the  Reynolds  stresses,  much  as  in 
the  previous  section.  This  is  a  natural  task  for  a  symbolic  manipulator  like  MACSYMA.  The  result  would 
enable  the  determination  of  all  unknown  coefficients  in  the  model  for  the  rapid  pressure  strain  term  (see 
Chapter  6);  we  are  attempting  to  carry  out  this  evaluation. 

4.7  Two-dimensional  turbulence 

RDT  of  two-dimensional  turbulence  is  useful  for  testing  the  range  of  performance  of  turbulence  mod¬ 
els.  Stanford  student  Laura  Pauley  carried  out  an  RDT  analysis  of  initially  axisymmetric  two-dimensional 
turbulence  for  three-dimensional  irrotational  strain  along  the  principal  axes  of  tl;,  with  62 2  =  -1/3.  Her 
results  are 

=  ^-fl  +  |ui  +  2a2  +  as]  +  —  (a^  +  a3)  +  aia3  +  2  a2(a[  +  a2  +  a3) 


11  3  5  2  11  2 


"*■  <*2(ai  +  as  +  2aia3)  +  25^  ^a!  + 


2-  ,  > 
3“2  +  “3/ 


+  0(a4) 


1  +  t(“‘  _  +  I^a'a3  “  a‘a^  +  °(a4 


(4.7.2) 


where  the  total  strain  in  the  i**  direction  is  =  exp(aj),  a2  =  a2  -  ao  and  ao  =  ai  +  a2  +  03.  Note  that  strain 
aligned  with  the  vorticity  does  not  affect  the  anisotropy,  and  that  changes  in  anisotropy  do  not  occur  until 
third  order.  It  would  be  instructive  and  useful  to  extend  this  analysis  to  more  general  2-D  cases  including 
rotation  and  strain  not  aligned  with  the  principal  axes  of  i>,;. 


5.  MODELING  SCALE  EVOLUTION  IN  HOMOGENEOUS  TURBULENCE 


5.1  Introduction 

This  chapter  and  the  next  are  devoted  to  one-point  models  of  homogeneous  turbulence.  Here  we  deal 
with  modeling  the  evolution  of  the  length  and  time  scales,  assuming  that  whatever  must  be  known  about 
the  tensor  character  of  the  turbulence  can  be  generated  by  an  anisotropy  model.  Anisotropy  modeling  is 
addressed  in  the  subsequent  chapter. 

The  turbulent  kinetic  energy  equation  provides  the  equation  for  the  turbulence  velocity  scale  q1.  For 
homogeneous  turbulence  (2.0.2)  becomes 

(42)  =  2(/>-e).  (5.1.1) 

If  the  other  model  scale  variable  is  et  this  equation  is  closed.  Alternatively,  if  one  choses  to  use  a  time  scale 
variable  r  instead,  a  model  relating  e  to  r  is  required.  There  are  many  clues  that  the  use  of  a  time  scale 
as  the  second  scale  parameter  would  offer  advantages  in  modeling  more  general  flows.  Where  it  becomes 
desirable  to  think  along  these  lines  we  will  use  the  large-eddy  time  scale  identified  in  Chapter  3, 

r  =  j-  (5.1.2)| 

Following  the  most  popular  current  trend,  we  shall  start  by  using  a  model  equation  for  £  as  our  second 
equation.  In  homogeneous  turbulence  (see  3.6.4)  e  =  i/w2,  and  so  the  iu2  equation  (2.8.12)  yields  the  e  (care 
must  be  taken  to  account  for  the  density  change  when  introducing  i/).  For  homogeneous  turbulence  with 
p  =  p(t)  and  n  =  constant  the  result  is 

e  =  2t/w1'ui'S1*  -  ^eS**  +  20,*;^'  +  2t'wl'ur's'J  -  2t/2w(„  u[,j  (5.1.3) 


is  the  anisotropic  strain  rate.  The  last  two  terms  provide  the  means  by  which  £  changes  in  isotropic  turbu¬ 
lence.  In  addition,  we  see  that  incompressible  strain,  isotropic  volume  change,  and  rotation  will  also  modify 
the  evolution  of  e.  We  shall  address  these  issues  separately. 

5.2  Decay  of  Isotropic  turbulence 

With  no  production  the  energy  equation  gives 

q7  =  -  2e  (5.2.1) 

Assuming  that  one  can  make  a  model  using  only  q2  and  e  as  variables,  the  form  of  the  e  equation  for  isotropic 
turbulence  can  be  deduced  by  dimensional  analysis, 

£2 

i=-C,  o—  (5.2.2) 


where  the  coefficient  C,o  can  depend  on  the  turbulence  Reynolds  number  (see  3.2.2)  Rt  =  qi/[i/e).  This  is 
the  form  used  by  all  models  of  this  type. 

Insight  is  obtained  by  recognizing  that  the  right  hand  side  of  (5.2.1)  comes  from  the  difference  of  the 
last  two  terms  in  (5.1.3).  The  first  of  these  is  the  turbulent  vortex  stretching  term,  which  is  related  to  the 
derivative  skewness  by  (3.10.9).  The  last  term  can  be  written  as 


2i/2<v',,  <u',y  =  2u2  [°°  k*E{k)dk 

Jo 


which  shows  that  it  is  dominated  by  the  smallest  scales  of  motion  and  hence  should  scale  on  the  Kolmogorov 
variables.  It  can  be  estimated  using  the  model  spectrum  of  Fig  3.9.4.  Using  these  two  estimates  for  isotropic 


turbulence,  put  in  terms  that  res  mble  the  model  equation 
and 


(5.2.1),  one  finds  that  both  terms  scale  as  vrRr  i 


10  \3or/  \q2' 


(52.4) 


Experiments  clearly  indicate  that  a  constant  coefficient  C,c  does  a  very  adequate  job  at  higu  Reynolds 
numbers,  which  means  that  the  difference  in  the  two  terms  within  the  brackets  in  (5.2.4)  must  decrease  as 
l/y/Rr-  Both  terms  are  very  large  and  they  are  nearly  in  balance  (an  estimate  of  the  skewness  can  be  made 
from  this  balance).  It  would  be  unwise  to  model  these  two  large  terms  separately  when  we  only  need  their 
difference,  and  for  this  reason  the  two  are  lumped  together  in  (5.2.2). 

The  value  of  C,0  can  be  determined  by  fitting  the  energy  decay  rate  for  isotropic  turbulence  to  that 
measured  experimentally,  and  this  is  what  most  modelers  have  done.  The  exact  solution  of  the  q2  and  e 
model  equations  is 

<j2  =  92(l-M/a)"  e  =  <r„(  1  +  t/a)  * 1 1  (5.2.5a, 6) 


a  = 


ffL 

2fn 


2 

(C. o  -  2) 


(5.2.5c,  d) 


The  subscripts  0  denote  initial  values.  The  best  experiments  suggest  n  should  be  in  the  range  1 . 1-1.3.  At  low 
Reynolds  numbers,  wl.^re  the  turbulence  is  in  its  final  period,  n  =  5/2  is  found  theoretically  and  confirmed 
experimentally. 

The  model  spectrum  (3,9.4)  cm  be  used  to  find  n  (Reynolds  1976)  by  assuming  that  the  spectrum  is 
permanent  below  fcp,  i.e.  that  the  low  wavenumber  spectrum  parameter  A  is  constant.  Expressing  i  in 
terms  of  q 2  and  A  using  (3.9.5)  and  (3.9.6),  then  using  this  in  (5.2.1)  to  find  the  q 2  history,  one  obtains 
(5.2.4a)  with  n  =  (2m  +  2)/(m+  3).  This  clearly  supportes  the  idea  that  the  low  wavenumber  part  of  the 
spectrum  affects  the  energy  decay  rate.  The  fc4  spectrum  (m  =  4)  gives  n  =  10/7,  which  is  really  too  high 
to  fit  the  best  experiments  very  well.  However,  the  fc2  spectrum,  with  m  =  2,  gives  n  —  6/5,  in  quite  good 
agreement  with  experiments.  In  a  finite- Fourier  series  representation,  the  assignment  of  the  same  energy  to 
each  low  wavenumber  Fourier  mode  would  make  E„  independent  of  k  and  hence  E[k)  vary  like  fc2,  and  so 
fc2  turbulence  can  be  thought  of  as  being  equipartitioncd  at  low  wavenumbers. 

With  n  =  6/5  as  suggested  by  both  the  experiments  and  the  fc2  spectrum,  C,0  =  11/3,  and  this  is  the 
value  that  we  prefer.  It  is  very  close  to  the  value  of  3.84  used  by  many  fc  -  e  modelers. 


5.3  Isotropic  compression 

For  isotropic  turbulence,  R<j  =  q7t>ij/ 3.  Denoting  5**  =  3r  (see  4.4.1),  and  assuming  isotropic  volume 
change  with  p  =  p(t),  the  energy  equation  (2.6.2)  reduces  to 


q2  =  — 2Tg2  —  2e.  (5.3.1) 


The  e  must  be  modified  to  account  for  the  change  in  volume.  The  exact  s  equation  (5.1.3)  suggests  that  this 
modification  might  be 


r  r 

e  =  -Ct0-~  -  el. 

<7 

For  very  large  T  the  solutions  to  the  above  equations  are 

(5.3.2) 

„2  __  2  — 2r< 

—  <7oc 

(5.3.3a) 

S  =  Co<-ri 

(5.3.36) 

The  energy  development  matches  RDT  (4.4.9).  If  we  assume  that  the  integral  scale  is  proportional  to  q3/e, 
the  large-eddy  length  scale,  then  according  to  (5.3.3)  the  length  scale  varies  as  exp(-2rt).  This  says  that 
expanding  the  flow  volume  will  reduce  the  length  scale,  which  should  be  disturbing  to  anyone  and  is  not  in 
agreement  with  RDT.  Nevertheless,  this  modification  of  the  e  equation  was  used  for  some  time  in  i.c.  engine 
modeling  before  the  problem  was  noted  (Reynolds  1980). 


The  RDT  analysis  suggests  instead  that  the  t  equation  for  this  problem  should  be 


(5.3.4) 


For  rapid  volume  change  this  produces 

e  =  c0c~*rt  (5.3.5) 

for  which  the  length  scale  varies  in  proportion  to  the  strain,  i.e.  as  exp(Tt). 

This  example  points  out  the  pitfalls  of  using  the  exact  equation  for  c  as  the  basis  for  its  model  equation. 
To  paraphrase  Saffman,  one  should  model  the  Physics  and  not  the  equations. 

5.4  Rotation 

Experiments  and  numerical  simulations  show  that  rotation  does  not  appreciably  alter  the  anisotropy  of 
isotropic  turbulence.  RDT  (section  4.6)  showed  that  rotation  does  not  affect  much  of  the  spectrum  at  all, 
but  does  tend  to  produce  a  slow  growth  in  the  energy  of  the  two-dimensional  component  of  the  turbulence 
aligned  with  the  axis  of  rotation.  The  simulations  (Bardina  et  al  1985)  reflect  this  growth  as  a  change  in 
the  integral  scales,  with  the  scale  in  the  direction  of  the  rotation  axis  becoming  longer  than  the  other  two  as 
time  passes.  Rotation  also  reduces  the  dissipation  rate,  apparently  by  inhibiting  the  energy  transfer  cascade. 

Most  turbulence  models  in  use  today  show  no  effect  of  pure  rotation  on  c,  a  weakness  that  has  been 
slow  to  receive  correction.  Bardina  found  that  his  large-eddy  simulations  and  Wiegland  and  Nagib's  (1978) 
experimental  data  could  both  be  predicted  extremely  well  using  a  simple  modification  of  the  c  equation, 


Hi! 

3  q2 


-  C.ijefl 


(5.4.1a) 


where  Q  is  the  rins  rotation  rate 

n  =  v/n.yfl.j.  (5.4.16) 

Bardina  found  that  Ctu  =  0. 15/ %/2  worked  well,  and  we  adopt  this  value. 

The  imposition  of  a  mean  strain-rate  provides  a  source  of  turbulent  kinetic  energy  through  the  turbulence 
production  term  (2.6.3).  'Ve  assume  that  the  anisotropy  part  of  the  turbulence  model  will  produce  f?,7 
values  given  q 2  and  e,  hence  P  need  not  be  modeled.  Thus,  no  modeling  for  the  q2  equation  is  required  for 
homogeneous  turbulence. 

The  associated  changes  in  the  dynamics  of  e  must  be  incorporated  in  the  e  model  equation.  To  date 
the  most  effective  means  for  doing  this  is  to  add  a  term  proportional  to  P, 


c2  Pc 

e  =  -C,0-z  -  C, nef)  +  C,p  — 

r  72 


(5.4.1) 


An  estimate  of  C,p  can  be  made  using  the  homogeneous  shear  flow  data  of  Tavoularis  and  Corrsin  (1981). 
Homogeneous  shear  flow  apparently  reaches  an  equilibrium  structure  in  which  the  Reynolds  stresses  all  scale 
with  the  turbulent  kinetic  energy.  The  energy  and  dissipation  rate  both  increase  with  time  in  a  manner  that 
keeps  the  turbulence  time  scale  very  nearly  constant  at  a  value  set  by  the  mean  shearing  rate  T  =  dUi/dz2. 
The  equation  for  r,  derived  from  (5.1.2)  using  (5.1.1)  and  (5.4.1),  is 


r  =  (C. o  -  2)  +  C. nflr  -  [C.p  -  2)  j.  (5.4.2) 

The  experiments  gave  T q2/c  =  12.7,  corresponding  to  fir  =  8.98,  and  P/c  =  1.8,  Using  C,n  =  11/3  and 
C,n  =  0.15/v/2,  a  constant  value  of  r  requires  C,p  =  3,45.  This  is  somewhat  higher  than  the  value  that 
Bardina  recommended,  which  was  based  on  his  large  eddy  simulations  of  strained  flows. 

Most  k  -  t  models  used  today  do  not  include  the  C, n  term.  For  plane  shear  flows  the  rotation  term 
and  the  production  terms  have  the  same  form,  and  when  these  terms  are  into  a  single  term  expressed  as  in 
the  form  of  the  production  term  the  resulting  combined  coefficient  based  on  the  above  coefficients  is  about 
3.0,  which  is  very  close  to  the  value  of  2.88  used  in  many  k  —  e  models. 
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S.5  Proposal  for  a  simple  k  —  r  model 

Compared  to  the  e  equation  (5.4.1),  the  r  equation  (5.4.2)  is  impressive  in  its  simplicity.  When  one 
examines  models  for  inhomogeneous  turbulence,  q7/e  frequently  appears,  suggesting  that  the  time  scale 
might  be  preferable  to  t  as  the  second  model  variable.  The  choice  should  be  based  on  the  ease  with  which 
the  model  extends  to  new  situations.  The  diffusion  terms  required  for  inhomogeneous  flows  are  particularly 
useful  in  evaluating  various  proposals. 

For  example,  consider  what  happens  to  the  terms  in  the  t  equation  near  a  solid  boundary.  The  t2  jq2 
term  goes  to  infinity,  but  the  Pefq2  term  goes  to  sero.  Consequently,  a  great  deal  of  effort  has  been  spent 
inventing  near-wall  patches  for  these  terms.  One  does  not  excape  these  simply  by  changing  variables,  unless 
a  slight  modification  is  made.  In  contrast,  Wilcox  (1986),  who  uses  a  reciprocal  time  scale  in  place  of  t, 
achieves  reasonable  near-wall  solutions,  even  in  the  viscous  region,  with  no  near-wall  modifications  of  his 
model  equation. 

Two-equation  models  have  been  criticised  because  the  length  scales  are  anisotropic  in  anisotropic  tur¬ 
bulence  but  the  model  assumes  isotropy  of  length  scale.  The  success  that  two-equation  models  enjoy  would 
seem  remarkable  in  the  light  of  this  objection.  But  suppose  it  is  really  fime  scale  information  that  is  carried 
by  e,  and  that  the  anisotropy  of  length  scales  is  reflected  by  anisotropy  of  R,,. 

Another  clue  is  provided  by  the  case  of  isotropic  volume  change,  for  which  the  r  equation  is 

r  =  (C,o  -  2)  -  (5.5.1) 


Note  the  appearance  of  the  strain  rate  term. 

It  is  suggested  that  it  might  be  better  to  replace  the  production  term  in  the  r  equation  by  a  term 
proportional  to  the  rms  strain  rate,  Some  additional  simplicity  of  form  is  obtained  by  using  the  kinetic 
energy  k  =  q2 / 2  and  redefining  the  time  scale  and  turbulent  Reynolds  number  by 

f  =  k/t  Rj  =  krfv  (5.5.2a,  4 ) 


The  model  equation  proposed  is 


f  =  Cm  +  C, oOf  -  Crs‘S‘f  - 

(5.5.3) 

Here  S’  is  the  rms  anisotropic  strain  rate 

S’  =  Js’  S’ 

V  'i 

(5.5.4a) 

determined  from  the  anisotropic  strain  rate  tensor 

S-,  =  Sii  -  \skk6t]. 

(5.5.44) 

Note  that  none  of  these  terms  is  ill-behaved  at  the  wall,  and  so  there  is  hope  that  the 
can  be  much  simpler.  The  constants  for  this  model,  evaluated  in  the  same  manner  as 
are 

C,  o  =  5/6  CVo  =  0.11  CrS-  =  0.69. 

near-wall  modifications 
those  in  the  c  equation, 

(5.5.5a,  6,  c) 

Exploration  of  this  idea  is  encouraged. 

6.  MODELING  ANISOTROPY  IN  HOMOGENEOUS  TURBULENCE 
6.1  Description  of  anisotropy 

The  scale  equations  developed  in  the  previous  chapter  are  closed  only  for  isotropic  turbulence  In 
general,  the  Reynolds  stress  tensor  must  also  be  determined  by  the  model.  The  Reynolds  stress  anisotropy 
tensor 

=  (6.U) 

9 

is  a  very  convenient  way  to  describe  the  deviations  from  isotropy.  This  chapter  deals  with  6, y  and  its 
modeling,  which  must  be  done  with  great  care  if  unrealistic  predictions  are  to  be  avoided. 

The  anisotropy  tensor  has  some  important  properties  that  need  to  be  kept  firmly  in  mind.  By  definition 
it  is  trace  free, 

fc„  =  0.  (6.1.2) 

It  is  often  convenient  to  think  of  bl}  in  its  principal  coordinates,  where  only  diagonal  elements  are  non-zero. 
By  (6.2. 1)  the  sum  of  these  principal  values  is  zero,  so  only  two  are  independent.  This  means  that  the 
anisotropy  can  be  characterized  by  two  independent  invariants, 


II  =  -btJbji/2 


III  = 


(6.1.3a,  b) 


If  the  turbulence  is  two-dimensional ,  meaning  that  one  (principal-axis)  velocity  component  is  always 
zero,  by  the  definition  (recall  Greek  indices  are  not  summed) 


ba„  =  -1/3 


if  Ra 


And,  if  all  of  the  energy  becomes  concentrated  in  one  component, 

6, ,,.  =  2/3  if  R„a=92.  (6.1.5) 

This  is  called  one-dimensional  turbulence.  Note  that  the  one  non-zero  velocity  component  could  be  a  function 
of  the  other  two  coordinates,  say  u',  (x2, 13,  J),  so  that  the  flow  would  resemble  a  honeycomb  of  opposing  jets. 

Thus,  the  possible  values  of  the  two  independent  principal  baa,  say  tjj  and  622,  must  lie  within  the 
triangle  on  Fig  6.1.1.  The  vertices  correspond  to  the  three  possible  states  of  one-dimensional  turbulence, 
and  the  sides  to  states  of  two-dimensional  turbulence.  The  isotropic  state  is  the  origin.  The  diagonal  lines, 
along  which  two  principal  components  are  the  same,  are  states  of  azisymmetric  turbulence. 


<-1/3, 2/3) 


<-1/3, -1/3) 


<2/3,-1/3> 


Figure  6.1.1  Range  of  possiole  principal  values  of  the  anisotropy  tensor 
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Note  that  one  can  either  move  along  an  axisymmetric  line  away  from  isotropy  to  a  two-dimensional  state 
(edge)  or  to  a  one-dimensional  state  (vertex).  These  limiting  cases  may  seem  extreme.  However,  turbulence 
near  a  wall  is  two-dimensional  (the  normal  component  vanishes),  and  turbulence  in  a  strongly  sheared  layer 
moves  remarkably  far  towards  one-dimensionality. 

In  homogeneous  turbulence,  the  move  towards  a  two-dimensional  state  is  made  by  straining  the  tur¬ 
bulence  in  one  direction  and  contracting  it  equally  in  the  other  two.  This  stretches  vortex  filaments  in  the 
direction  of  positive  strain,  aligning  these  filaments  with  the  flow  and  thereby  reducing  the  fluctuations  in  the 
direction  of  positive  strain.  This  is  what  happens  to  turbulence  when  it  is  passed  through  an  axisymmetric 
contraction. 

The  move  towards  a  one-dimensional  state  is  achieved  by  straining  the  flow  equally  in  two  orthogonal 
directions,  and  contracting  it  in  the  third,  as  one  could  do  in  an  axisymmetric  diffuser  (using  boundary 
layer  suction  to  prevent  separation).  The  vortex  cores  are  stretched  out  to  form  sheets  (pancakes)  and  the 
limiting  one-dimensional  case  corresponds  to  a  honeycomb  of  two-dimensional  vorticity.  We  will  call  this 
type  of  deformation  axisymmetric  expansion. 

An  equivalent  and  less  specific  way  to  characterize  the  anisotropy  is  through  the  anisotropy  invariant  map 
introduced  by  Lumley.  For  axisymmetric  turbulence  we  write  the  anisotropy  tensor  in  principal  coordinates 

0  0  \ 

a  0  .  (6.1.C) 

0  -2  a  J 

III  = -2a3.  (6.1.7) 


Along  lines  where  a  <  0  so  that  the  component  along  the  axis  is  more  energetic  than  the  other  two  (axisym- 
metric  expansion), 


III  =  +2 


(6.1.7a) 


while  if  a  >  0  so  that  the  axis  component  is  less  energetic  (axisymmetric  contraction) 


III  —  —2 


3/2 


(6. 1.7a,  b) 


The  two-dimensional  boundaries  can  be  studied  in  principal  coordinates,  writing 


0 

Oil 

c 

0 


Then 


so  that  it  for  two-dimensional  turbulence 


III  = 


a2  —  1 
108 


G  =  ^  +  11  +  3111  =  0. 


(6.1.8) 


(6.1.9a,  b) 


(6.1.10) 


Using  these  results,  the  range  of  possible  turbulence  states  is  shown  in  the  invariant  map  of  Fig.  (6.1.2) 
The  origin  is  the  isotropic  state,  the  upper  boundary  is  the  locus  of  two-dimensional  states,  the  two  sides 
are  the  two  types  of  axisymmetric  states,  and  the  upper  vertex  is  the  one-dimensional  state.  The  anisotropy 
invariant  map  is  a  very  useful  way  to  characterize  the  state  of  turbulence  in  modeling,  simulations,  and 
experiments. 

Two  tensors  that  can  be  formed  from  the  anisotropy  tensor  are  its  square, 


i>?j  =  bikbkj 


(6.1.11) 


l-4*i 


and  its  cube, 

6^  =  6,„6„m6m/-  (6.112) 

The  Cayley- Hamilton  theorem  of  linear  algebra  says  that  a  matrix  satisfies  its  own  characteristic  equation, 
which  in  this  instance  means  that 

b*j  +  Ittjy  -  IlWjy  =  0  (6.1.13a) 

or  alternatively 

»?/  =  \blkb„  +  l-b3kkS„.  (6.1.136) 

Hence,  6fy,  and  all  higher  powers  of  the  tensor,  are  linearly  dependent  on  the  lower  powers  and  hence  do 
not  contain  new  tensorial  structure  beyond  that  in  6?  ,  6,y,  and  6,y.  As  we  shall  see,  this  is  very  important 
in  turbulence  modeling.  Readers  not  familiar  with  this  important  theorem  may  find  it  instructive  to  verify 
(6.3.13b)  by  writing  6<y  in  its  principal  coordinates,  carrying  out  the  products  using  the  trace-free  condition 
to  express  one  of  the  principal  values  in  terms  of  the  other  wo. 


(2/27,1/3) 


Figure  6.1.2  Anisotropy  invariant  map 


6.2  Evolution  equation  for  the  anisotropy  tensor 

Using  the  evolution  equation  for  R< ,  (2.7.1)  and  the  definition  of  6,-y,  the  equation  for  evolution  of  6,, 
in  homogeneous  turbulence  with  p  =  p(t)  can  be  written  as 


bij  —  {bikSkj-  +  bji,Ski  ,bnmSnmSjJ)  +  26n(n5nm6ty 

+  (6i*H*y  +  b]k  f}*,)  -+■  —  [T,y  —  (Dij  —  “  Dkk6tl )  ]  -t-  2— 6,-y 
g  a  g 


(6.2.1) 


Here  S'y  is  the  anisotropic  strain  rate  tensor  defined  by  (5.5.4).  Note  that  only  anisotropic  strain  produces 
Reynolds  stress  anisotropy,  and  that  the  right  hand  side  is  properly  trace-free.  The  terms  containing  the 
mean  rotation  tensor  ft,/  represent  a  kinematic  rotation  of  the  anisotropy  tensor.  When  used  in  conjunction 
with  q1  and  e  equations,  models  for  the  pressure-strain  term  T.j  and  the  anisotropy  of  the  dissipation  tensor 
Dij  must  be  provided. 

There  is  a  class  of  turbulence  models  called  algebraic  two-equation  models  in  which  iv  ><•  assumed  that 
the  turbulence  structure  has  reached  an  equilibrium  state  determined  by  a  balance  of  the  terms  on  the  right 
hand  side  of  (6.2.1).  For  example,  the  standard  k  -  e  model  uses  an  algebraic  equation  equivalent  to 

bi,  =  (6.2.2) 

with  =  0.09.  A  problem  should  be  immediately  apparent.  The  sudden  imposition  of  a  strong  strain 
could  easily  produce  6,/  states  lying  outside  of  the  anisotropy  invariant  map.  This  is  a  very  serious  potential 
problem  when  such  models  are  applied  in  new  flows. 

Another  weakenss  of  this  model  is  that  it  assumes  that  the  principal  axes  of  stress  and  strain-rate  are 
aligned.  This  is  not  true  in  the  most  important  engineering  flow,  namely  shear  flow.  However,  the  constant 
has  been  set  to  give  the  right  anisotropy  of  the  shearing  stress.  For  example,  in  the  homogeneous  shear 
flow  experiment  of  Tavoularis  and  Corrsin  discussed  in  section  (5.4),  bt 2  =  —0.149  is  predicted  by  (6.2.2),  in 
excellent  agreement  with  the  measurements.  However,  the  model  predicts  fcjj  =  0,  whereas  the  experiments 
show  6U  =  0.196,  so  the  normal  stresses  are  badly  in  error.  However,  they  do  not  play  a  significant  role  in 
determining  the  mean  velocity  field,  and  so  this  error  usually  of  little  consequence. 

Algebraic  models  assume  that  the  turbulence  structure  responds  instantly  to  changes  in  the  imposed 
mean  strain.  This  is  reasonable  for  computing  the  slow  evolution  of  mean  fields,  but  not  satisfactory  if  the 
strain  rates  are  large,  i.e.  if  S'q^/c  0,  where  S'  is  the  rms  anisotropic  strain  rate.  And  algebraic  models 
predict  instant  restoration  of  isotropy  after  the  removal  of  an  applied  mean  Btrain-rate.  Hence,  if  one  wants 
to  have  realistic  predictions  of  the  Reynolds  stresses  in  these  cases,  a  model  of  the  6,/  evolution  equation  be 
solved  in  parallel  with  the  q2  and  e  equations. 

In  the  balance  of  this  chapter  we  review  the  formal  methods  that  have  been  applied  in  attempts  to 
develop  rational  models  to  close  the  6,,  evolution  equation.  Then,  at  the  end  we  will  present  a  much  simpler 
model  that  achieves  some  of  the  objectives  of  the  more  complicated  models  at  much  less  expense.  This  new 
model  might  be  useful  for  engineering  analysis. 

6.3  Decomposition  of  the  pressure-strain  term 

The  Poisson  equation  for  the  fluctuation  pressure  (2.6.1)  has  two  terms  on  the  right  that  act  as  sources 
for  pressure  fluctuations.  The  source  involving  the  mean  velocity  gradients  will  change  instantly  when  the 
gradients  change,  resulting  in  an  instant  change  in  the  fluctuating  pressure  field  and  hence  an  instant  change 
in  the  pressure-strain  term  T,/.  The  source  involving  only  the  turbulence  will  change  only  as  the  turbulence 
adjustes  to  its  new  consitions.  This  suggests  that  the  pressure  fluctuations  be  split  into  rapid  and  slots  parts, 

p'=p*‘*  +  p*?)  (6.3.1a) 

where  the  rapid  term  is  the  solution  of 

-P,l|,i.  =  ~2u'  „  f/,,y  (6.3.16) 

P 


and  the  slow  term  is  the  solution  of 


1 

~P 

P 


(2’,.i  =  -2 u'„  uj „•  +2u'.„  uj, j. 


(6.3.1c) 


The  resulting  contributions  to  the  pressure-strain  term  (2.7.5)  will  be  denoted  by  T-p  and  \  respectively. 

Eqn.  (3.6.1b)  is  linear  and  has  constant  coefficients  in  homogeneous  turbulence,  and  so  can  be  solved 
by  Fourier  methods.  We  follow  the  approach  of  Chapter  3,  and  write 


\ 
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P«‘»  =  x>(k)e'“‘ 


(6.3.2a, t) 


The  time-dependence  of  the  coefficients  is  not  explicitly  expressed  because  we  are  solving  the  Poisson  equation 
at  one  instant  of  time.  The  solution  of  (6.3.1b)  is  then 

-p(k)  =  -2 Ur„  ^fiy(k).  (6.3.3) 

p  K 

Multiplying  the  pressure  fluctuation  series  by  the  the  velocity  gradient  series,  using  the  conjugate 
symmetry  properties  of  the  Fourier  inodes,  averaging  over  the  box  of  Fig.  3.3.1,  then  taking  the  limit 
as  done  in  section  3.3,  one  finds 

=  2 UPi]  Mt]vq  (0.3.4 ) 

P 


M,„„I  =  J  ^pp£,y(k)d3k. 


The  rapid  pressure-strain  term  is  the  sum  of  two  such  terms, 


Tiq  —  21/,,,/  +  A/,,/,„).  (C.3.C) 

Modeling  of  the  rapid  pressure  strain  term  therefore  becomes  a  task  of  modeling  A /,,,„,,  which  we  address 
in  the  next  section. 


6.4  Modeling  the  MiJvq  tensor 

The  Mi/,,,  tensor  has  been  modeled  in  various  ways,  all  relatively  simple,  usually  with  one  constant 
being  adjusted  to  fit  data  for  the  predictions  of  a  selection  of  flows.  Here  we  introduce  a  very  different 
approach;  we  arge  that  the  anisotropy  model,  when  applied  in  circumstances  for  which  rapid  distortion 
theory  would  apply,  should  give  results  consistent  with  RDT.  The  RDT  form  of  the  model  equation  includes 
only  the  rapid  pressure  strain  term,  the  production  term,  and  the  mean  rotation  term  in  (6.2.1),  exactly 
the  same  terms  used  in  RDT  theory.  The  solution  above  for  the  rapid  pressure  field  is  exactly  the  same  as 
used  in  RDT.  Therefore,  in  principle  it  should  be  possible  to  determine  all  of  the  coefficients  in  the  rapid 
pressure  strain  model  (i.e.  in  A/,/,,,)  so  as  to  make  the  anisotropy  predicted  by  the  model  equation  under 
RDT  approximations  exactly  the  same  as  that  predicted  by  RDT  theory,  for  an  arbitrary  rapid  deformation. 

Following  Shih  and  Lumley  (1985),  we  begin  by  writing  the  general  expression  for  a  tensor  A/,/,,,,  = 
Afy, ,'Jq2  that  is  assumed  to  be  a  function  of  the  tensor  6,,,  with  the  symmetries  in  indices  required  by  the 
definition.  This  is 


A/,/, i,  Ci  6tj  <5,,,/  d-  C2  (bii>byq  ~b  /,,)  -f-  C3 +  I'l b[itj b , j  -f"  f-s(tfi,,6/,  -f  6t  ijbjfi  “b  bj ,,6,,,  d-  bj ,, /*,,,) 

+Co6ijbf,q  +  CjSMbfy  +  C8  (6,,,bj,t  +  6i,,b2v  +  6,rbfq  +  <5//,/,)  +  Cgh./A,,,,  +  Cio{bi,,bjq  +  h,v6/,,) 

+Cnbijblq  Cxibpqbfj  +  C13(6ip6/,  +  b,qb2p  +  bjpbfq  +  b,qb2p)  +  C146*  bpq  +  Cis(b?pbfq  +  b2qb2lp).  (6.4.1) 

Because  of  the  Cayley-Hamilton  theorem,  higher  powers  of  6,/  are  not  required.  The  coefficients  C ,  -  C15 
may  be  function  of  the  invariants  II  and  III,  and  of  other  scalars,  such  as  Rt- 

The  continuity  (2.3.4)  equation  requires  A/,/,,  =  0.  When  this  condition  is  applied  to  (6.4.1),  an 
equation  containing  6jq,  bjq,  and  b2q  is  obtained.  Since  these  are  independent  tensors,  the  coefficient  of  each 
must  vanish.  This  produces  three  equations, 

Ci  +  4C2  +  b2kCg  +  -bkk[Cn  +  Ci  2  +  2C13)  =  0  (6.4.2a) 

C3  +  Cq  +  5C5  +  -&fc*(C|i  +  C12  +  4  C,3)  +  -bkk(Cu  +  C15)  =  0 


(6.4.26) 
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Cg  +  Ct  +  SC/  +  Co  +  Cjo  +  — +  SCit)  —  0.  (6,4.2c) 

FVom  iti  definition,  Mi)pq  =  Jfcy,  so  Mitfyq  =  by  +  Sy/3.  In  the  «ame  way,  this  condition  gives  three 
additional  constraints, 

3Ct  +  2CS  +  blkCt  +  -b3kkCl3  =  -  (6.4.3a) 

3Ct  +  4  Co  +  bkk(Ci  i  +  2C13)  +  —fifcfcC'js  =  1  (6.4.36) 

3C7  +  4Cg  +  2Cjo  +  fcfcfc  (<?i  4  +  Ci6)  =  0-  (6.4.3c) 

These  six  conditions  reduce  the  number  of  undetermined  coefficients  to  nine,  and  give 

Ci  =  h + 4*‘  (sCs  +  \Ci  +  T&Ci0) +  tkk  (ic“ +  iCn  ~  !Cl3) 

*'(i**)3  +  s^li)  (6.4.4a) 

°3  =  +  i“(~!C8  -  T5C°  ~  Ec'°)  +  bkk{~Toc "  ~  BCl2  "  Is0'3 

+^{-EClt~ToCli)  {&AAb) 

cs  =  -  YCi  +  ("s^11  ~  2Cn  ~  J^13)  +  4**  (“i^14  "  gi^16)  (6.4.4c) 

=  5  -  ic* + tkk  {ACn  ~  lCi3) + iik  HCis)  (6AAd) 

C0  =  - yCs  -  Co  -  jCjo  +  b\k  ^--Cu  -  -C1S^  (6.4.4e) 

Cl  =  - \Co  -  jC10  +  blk  (-jCh  -  jCls)  (6.4.4/) 

Once  the  coefficients  are  evaluated,  T*.1*  can  be  determined.  The  result,  written  in  terms  of  the 
anisotropic  strain  rate  tensor  (5.5.4)  and  the  rotation  tensor,  is 

T-l>  (  2  \ 

=  2 (C,  +  Cj)5,*y  +  (C,  +  c«  +  2C5)  ys‘kbkl  +  S‘kbki  -  -s;mbmns„  J 

+(C0  +  c7  +  2 c8)(sx  +  s;^  -  55^6^6,,)  +  2(cs  4  cw)s;q(b,qbJP  -  +  2C10s;,6,p6,y 

+  (Cii  +  Cia  +  2C13)Sp*?  ^biqbp)-  +  6y,6^  -  -6^n6^J  +  2C13S^]bqpb,]  +  2Ci3Spqbqv  ^6,*-  -  -62„6,,^ 

+2(Cj«  +  CJ5)S^  (i2p62,  -  j62p62,6i;)  +  2C165p*,6^,(6?y  -  ji*.^) 

+(C3  -  C4)(flt,6k,  +  n*,6*)  +  (Ce  -  C7)(nkl6j|,  +  nky6?,)  +  (C„  -  Cn)np,(6„6^  +  4y,£).  (6.4.5) 

Realizability  has  been  of  much  concern  in  modeling  the  pressure-strain  term  and  other  terms  in  the 
bij  equation.  The  principal  values  baa  can  not  be  less  than  -1/3,  and  any  model  that  would  carry  a 
principal  value  below  this  amount  (i.e.  outside  the  bounds  of  the  invariant  map)  then  produces  unrealizable 
turbulence  (nonsense).  Truncated  approximations  to  the  series  above  have  this  danger,  although  the  model 


just  described,  with  the  infinite  set  of  coefficients,  would  be  realisable  because  RDT  solutions  are  realizable. 
In  order  to  guarantee  realizabilty  one  can  enforce  certain  conditions.  There  are  various  ways  to  develop 
these  conditions.  Shih  and  Lumley  (198S)  get  them  by  requiring  that  the  Taa  terms  must  not  drive  baa  out 
of  bounds.  This  requires 

Up,j  =  0  when  =  -1/3  (6.4.0) 

which  produces  three  additional  constraints, 

+  j(c’3+<?4+2C5)  +  -(Co+C7+2Cg  +  C,#+C’,o)-— (C,u+C,i2+2C13)  +  — (C,4  +  Ci5)  =  0  (6.4.7a) 


Cf,  —  ~<?10  gCl3  —  0 

C8  -  -C13  +  -C15  =  0, 


(6.4.76) 

(6.4.7c) 


We  believe  it  is  preferable  to  impose  the  rea.iz ability  conditions  directly  on  the  modeled  tensor  .  When 
the  velocity  component  u'a  is  everywhere  zero  then  u„  =  0  and  consequently  Majpq  =  0.  In  the  principal 
coordinates  of  6,y  this  requires 


Mim-0  ~  Mn 


M  1212  =  0  when  6u=— 1/3 


(6.4.8a,  6,  c) 


Using  the  fact  that  6jJk  =  -1/9  +  bkk/2  on  the  two-dimensional  line,  (6.4.8a-c)  give 


Ci  +  2C?  -  ~(C3  +  C4  +  4C5)  +  -(C6  +  C7  +  4Cs  +  Co  +  2Cjo)  -  —  (Cn  +  ^12  +4C|3)  +  —  (Ci«  +  2Cis)  —  0 

(6.4.9a) 

-C3  +  -(Co  -  Cg)  -  — (Cji  -  Cn)  +  —  Cm  =  0  (6.4.96) 


5Cs  +  i(C8_Clo)  +  Ii2Cl5  =  0 


(6.4.9c) 


C2  “  \Ci  +  (]S  +  \bkk)c*  -  T»Cl° +  (sr  ■  i6**)Cl3  +  (“152 +  Tshkk)Clb  =  °-  (6-4  M) 

When  equations  (6.4.4)  are  used  to  express  the  lower  coefficients,  (6.4.9d)  is  -1/2  (6.4.9a),  and  so  only  three 
independent  conditions  are  obtained, 

IS  +  5C* +  ("5  +  C<i  +  £**fcC* +  (J +  ihkk)c,° +  +  Tobkk) [Cn  +  c,2) 

+  (_45  +  5i*fc)C'13  +  3g(^*t)2^"14  +  +  •jq(***)*^ C'is  ~  0  (6.4.10a) 

~l  "  TC6  ”  Tl°*  ~  lC°  ~  ~  +  ^*)C“  +  (F  '  lb“)C" 

~lb2kkCi3  +  {rr  ib2kk)Cit+ { Wi~iibkk)Cib  =  o  (6-410t) 


ECs  +  T8(c°  _  Cio)  +  lfeCl6  =  0 


(6.4.10c) 


When  these  are  satisfied,  the  Shih-Lumley  conditions  will  also  be  satisfied.  It  is  important  to  realize  that 
these  realizability  constraints  apply  only  when  the  turbulence  is  two-dimensional,  i.e.  only  on  the  line  G  =  0 
that  forms  the  top  boundary  of  the  invariant  map. 

The  equations  above  suggest  that  the  coefficients  will  depend  on  the  invariants  and  not  simply  be 
constants.  We  might  expand  each  coefficient  as  a  power  series  in  the  invariants, 

c„  =  C<°»  +  b\kc™  +  b3kkcw  +  (&)aci4»  +  blkblkc^  +  (b3kkfcle>  + ... 


(6.4.11) 
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We  see  that  the  first  approximations  to  the  isotropic  coefficients  Ci  and  Cj  are  already  known,  and  the  first 
approximations  to  the  linear  coefficients  C3-C3  a re  determined  by  the  first  approximation  to  CE. 

Most  turbulence  models  presently  in  use  include  only  the  terms  in  through  C&  (the  linear  terms), 

employing  constant  values  for  the  coefficients.  But  with  Ce-C iE  =  0  no  single  value  of  CE  can  satisfy  all 
three  realisability  conditions  (6.4.10),  so  these  linear  models  do  not  satisfy  realizability. 

The  simplest  set  of  coefficients  satisfying  realizability  is  obtained  by  truncating  A4,JP,  to  0(b2)  and 
assuming  all  coefficients  are  constants.  The  truncation  gives 

Cu  =  Ci 2  =  C\3  =  C\ 4  =  C\ 5  =  0.  (6.4.12a) 

Prom  (6. 4.4a, b),  the  coefficients  Ci  and  C2  will  be  constants  only  if  C8  =  0  and  Cm  =  -3C9.  Then,  the 
realizability  conditions  give 

Ci  =2/15  C2  =  - 1/30  C3  =  1/30  C4  =  7/15 

C6  =  -l/10  C7  =  2/10  C9  =  1/10  C10  = -3/10  (6.4.126) 

These  are  the  coefficients  determined  in  a  slightly  different  manner  by  Shih  and  Lumley  (1986).  Under  the 
rapid  distortion  approximations,  the  time-series  solution  of  the  model  equations  resulting  from  (6.4.12)  match 
RDT  of  isotropic  turbulence  only  to  0(1).  The  model  also  predicts  that  anisotropic  turbulence  subjected  to 
pure  rotation  would  undergo  anisotropy  changes,  in  excess  of  those  caused  by  the  kinematic  rotation  terms, 
of  O(t),  whereas  RDT  indicates  that  this  excess  change  must  be  an  even  power  series  in  t  (see  section  4.3) 
and  hence  should  not  appear  until  0(t2).  It  would  seem  desirable  to  obtain  a  better  match  to  RDT. 

Under  rapid  pure  rotation  of  anisotropic  turbulence,  (6.4.5)  will  produce  an  O(t)  change  in  6i;  in  excess 
of  that  produced  by  the  kinematic  rotation  terms  unless  (C^01  -  <?{°!)  =  0.  This  condition  gives 

Cf1  =  -2/7  C<°!  =  C<0)  =  5/7.  (6.4.13) 

With  these  values,  the  RDT-equivalent  model  predicitions  also  agrees  with  RDT  to  O(t)  for  all  irrotational 
strains  (Reynolds  1983).  Le  Penven  and  Gence  (1983)  carried  the  analysis  to  one  additional  order  in  t  for  the 
case  of  irrotational  strain  at  a  constant  strain  rate,  and  found  that  the  coefficients  could  indeed  be  matched 
to  0(t2).  Hence,  it  seems  clear  that  (6.4.13)  gives  the  rational  choices  for  the  first  approximations  to  the 
linear  coefficients.  However,  with  CE  =  -2/7  the  realizability  conditions  can  not  be  satisfied  by  a  truncation 
of  Mi]pq  to  0(b2),  and  one  must  include  higher-order  terms  to  effect  realizability. 

It  seems  clear  that  continued  matching  with  RDT  would  determine  all  of  the  coefficients,  and  since 
RDT  predicts  realizable  turbulence  the  resulting  model  would  guarantee  realizability.  The  RDT  required  for 
a  complete  matching  must  be  sufficiently  genera)  to  allow  all  coefficients  to  be  det’rmined.  The  arbitrary 
irrotational  strain  analysis  given  in  section  (4.5)  is  not  sufficient  because  there  the  principal  axes  of  S,y  were 
fixed  and  hence  the  principal  axes  of  S,j  and  6,;  always  remained  aligned.  An  RDT  for  of  isotropic  turbulence 
with  arbitrary  initial  rotation  rate  and  arbitrary  strain  rate  history  is  required  (see  section  4.6).  It  should  be 
possible  to  select  the  constants  in  the  coefficient  expansions  (4.6.10)  to  match  RDT  to  any  arbitrary  order  in 
a  time-series  solution  of  the  RDT-approximate  model  equations,  and  then  to  use  the  realizability  conditions 
to  truncate  the  expansions,  maintaining  full  realizability.  Thus,  in  principal  the  rapid  pressure  strain  model 
should  be  determined  completely  by  RDT  analysis,  with  no  adjustable  constants  matched  to  experiments. 
We  are  attempting  to  complete  this  task. 

Another  approach  that  may  be  fruitful  is  to  use  RDT  for  initially  axisymmetric  two-dimensional  tur¬ 
bulence,  in  conjunction  with  the  realizability  constraints,  to  develop  expressions  for  the  coefficients  that 
must  hold  along  the  two-dimensional  line  G  =  0.  The  results  of  section  4.7  should  be  useful  in  this  regard. 
These  coefficients  might  then  be  expanded  in  power  series  in  G  in  order  to  determine  appropriate  values 
for  three-dimensional  turbulence,  perhaps  by  matching  to  RDT.  Many  interesting  analyses  of  this  nature 
remain  to  be  done  in  turbulence  modeling. 

6.5  Modeling  the  slow  terms 

The  negative  of  the  slow  pressure-strain  term  and  the  dissipation  ansiotropy  term  are  modeled  together 
in  (6.2.1)  as 

-  (D„  -  Dkk6,,/3)  =  (6.5.1) 
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Assuming  that  it  xs  possible  to  model  in  terms  of  i.y,  a  premise  that  is  not  supported  very  well  by  direct 
numerical  simulations,  the  most  general  form  must  be 


fa,  =  (a  +  2)fc,j  +  £(6*  +  —  IIA.j ) 


(6.5.2) 


where  the  coefficients  a  and  f)  could  be  functions  of  the  invariants  II  and  III  and  possibly  of  other  scalars,  such 
as  Rr-  Under  these  assumptions,  one  can  in  principle  evaluate  the  coefficients  by  reference  to  experiments 
and  simulations  on  the  return  to  isotropy  following  removal  of  mean  strain  rate.  In  this  case  (6.2.1)  reduces 
to 


i><j  =  — -  2b<i)  =  +  Wi \  ~  T^k4.;)]- 


(6.5.3) 


If  the  anisotropy  is  weak,  a  controls  the  return  and  must  be  positive  if  there  is  to  be  a  return. 

Using  (6.5.3),  the  evolution  of  the  state  point  on  the  invariant  map  is  described  by  the  two  equations 


dn 

dt 


-4(2oII-  3/9III) 

9 


(6.5.4) 


dill 

dt 


-4(3aIII+^II2) 

9  3 


so  that  the  trajectory  on  the  map  is  described  by 


dll  _  2aII-  3/3III 
dill  “  3aIII  +  |/?IIJ 


(6.5.5) 


(6.5.6) 


Therefore,  if  the  underlying  premise  of  the  model  is  correct,  the  trajectories  must  be  unique  and  the  ratio 
*7  (II,  III)  =  a//J  can  be  determined  by  the  local  trajectory. 

There  have  not  been  many  experiments  on  the  return  to  isotropy.  Those  that  do  exist  often  show  very 
strange  behavior.  Direct  numerical  simulations  of  Lee  and  Reynolds  (1985)  using  the  Rogallo  code  in  a 
1283  mesh  attempted  to  address  the  e  questions  in  the  hope  of  evaluating  the  parameters.  Turbulence  that 
had  been  strained  by  axisymmetric  contraction  relaxed  smoothly  to  isotropy  along  the  axisymmetric  line  as 
expected.  But  turbulence  that  had  been  strained  by  axisymmetric  expansion  showed  very  strange  behavior, 
in  some  cases  moving  further  away  from  isotropy  before  starting  the  return.  Turbulence  strained  by  complex 
combinations  that  produced  states  near  the  middle  of  the  anisotropy  map  did  not  show  convincingly  unique 
trajectories.  A  sample  of  the  trajectories  following  removal  of  plane  strain  are  shown  in  Fig.  6.5.1.  The 
points  to  the  left  have  been  strained  most  rapidly,  and  the  initial  states  are  preducted  very  well  by  RDT. 
The  lowermost  points  are  in  general  agreement  with  the  one  experiment  on  the  relaxation  from  plane  strain 
by  Tucker  and  A.  Reynolds  (1968).  Note  that  one  point  begins  its  “return”  by  going  substantially  far  in  the 
wrong  direction.  It  seems  impossible  to  incorporate  this  wierd  behavior  within  the  structure  of  (6.5.2). 


f 
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Figure  6.5.1  TYajectories  of  the  return  to  isotropy  from  plane  Btrain  (simulations) 

The  simulations  cast  doubt  on  the  basic  idea  of  modeling  these  terms  using  only  the  tensor.  Bu  . 
the  simulations  did  show  that  the  return  of  the  small-scales  to  isotropy,  as  reflected  by  the  anisotropy  in  the 
vorticity  and  disssipation  tensors,  was  quite  well  behaved  and  easily  modeled.  This  suggests  some  directions 
for  future  modeling  research. 

These  simulations,  as  well  of  those  of  Rogallo  for  homogeneous  shear  flow,  suggest  very  strongly  that 

4>ij  — >  26;y  as  |6|  — •  0.  (6.5.7). 

This  means  that  there  should  be  no  linear  return  to  isotropy.  Careful  examination  of  the  very  nearly  isotropic 
data  of  Comte-Bellot  and  Corrsin  (1966)  seems  to  support  this  behavior. 

Choi  (1983)  perfomed  experiments  on  the  return  to  isotropy  from  the  right  side  of  the  invariant  map, 
and  did  seem  to  observe  more  consistent  behavior.  A  fit  to  his  data  developed  by  the  Cornell  group  and 
reported  by  Shih  and  Lumley  (1986)  is 

a  =  12.44(9G)J(1  -  9G)3'*  0  =  0.  (6.5.8) 

The  G  factors  provide  a  sort  of  realizability,  and  there  is  no  linear  return  to  isotropy. 


1 

i 
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A  criticism  that  might  be  raised  about  this  model  is  that  it  does  not  allow  two-dimensional  turbulence  to 
remain  two-dimensiona',  relaxing  to  an  axisymmetric  state.  It  is  possible  to  construct  a  model  that  does  by 
using  the  realizability  c  ndition,  When  u'a  =  0  everywhere,  then  D„a  =  0,  Ta„  —  0,  and  hence  </„„  =  26u„, 
which  will  sustain  1/3.  Thus,  the  readability  condition  gives 


a  =  3/3 


when  G  =  0. 


(6.5.9) 


Since  this  constraint  only  need  be  true  for  G  =  0,  we  can  add  functions  of  G  without  destroying  realizability. 
A  linear  term  suffices,  with  its  coefficient  chosen  to  remove  the  linear  return  to  isotropy  when  (7=1/9  and 
II  =  0  and  to  make  /?  vanish  for  small  anisotropy, 


a  — 


Q  +  2II-3c)ft 


(6.5.10a) 


/3  =  ft(l  -9(7). 


(6.5.106) 


The  model  is  then 


^  j  +  211  -  3(7^oi’ij  +  Ai(  1  ~  9(7) 


(6.5.11) 


With  this  model,  for  nearly  isotropic  turbulence  (6.5.4)  becomes 


a  *  -ft  II 


(6.5.12) 


while  for  small  anisotropy  (6.5.8)  gives 


a  «  12.44(— 9II)3/4. 


(6.5.13) 


Matching  at  —II  =  0.05  suggests  ft  *»  10.  This  modified  model  satisfies  realizability,  restores  axisymmetry 
in  two-dimensional  turbulence,  displays  no  linear  return  to  isotroy,  and  gives  return  rates  of  the  right  order 
of  magnitude. 

However,  one  might  suspect  that  the  slightest  little  three-dimensionality  would  explode  the  turbulence 
into  a  three-dimensional  field,  so  perhaps  it  is  unreasonable  to  insist  on  maintaining  two-dimensionality  in 
the  model.  Undecided  issues  likes  this  provide  fruitful  grounds  for  new  research,  and  we  are  now  exploring 
questions  like  these  using  direct  turbulence  simulation. 
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8.6  A  simple  anisotropy  model 


The  gap  b  tween  the  eddy-viscosity  models  used  in  the  simplest  k-e  models  and  those  discussed  above 
is  immense.  There  is  a  need  for  a  much  simpler  model  that  would  protect  engineering  calculations  from 
the  dangers  of  unrealizable  turbulence,  provide  some  indications  of  the  trends  in  anisotropy  for  unusual 
flow  situations,  and  handle  dynamic  changes  on  roughly  the  right  time  scale,  but  without  such  calculational 
complexity.  The  beginnings  of  such  an  idea  are  -  jented  here. 

We  start  with  the  idea  that  a  large  positive  strain  rate  in  one  direction  tends  to  stretch  vortex  filaments 
in  that  direction,  aligning  them  with  the  flow,  thereby  intensifying  the  perpendicular  fluctuation  components 
and  reducing  those  along  the  axis.  In  the  limit  of  very  strong  strain  rate,  the  energy  in  the  axial  fluctuations 
axial  will  approach  zero.  The  anisotropy  model  must  prevent  negative  values.  And,  we  know  that  only  the 
it  anisotropic  component  of  strain  produces  anisotropy  in  the  turbulence.  A  simple  algebraic  model  with 
this  character  is 


>>.,= 


SIil 


Ao  4-  AsS*t 


(6.6.1) 


In  order  for  realizability  to  be  maintained,  baa  should  approach  -1/3  as  S'a  — >  oo,  for  any  combinations  of 
other  S'- .  This  requires  that  the  coefficient  As  depend  on  the  type  of  strain. 

In  the  principal  coordinates  of  S'  ,  we  take  S as  having  a  large  positive  value  I',  and  write  the 
strain-rate  tensor  as 

/!  0  0  \ 

S*=r  0  -ift  0  .  (6.6.2) 

\o  o  -^r) 


Note  that  a  —  0  gives  axisymmetric  contraction,  a  —  1  gives  plane  strain,  and  a  =  3  gives  axisymmetric 
expansion.  Then 


S' 


=  rv^ 


1  a.  t1  +  a)2  ,  I1  ~  a)2  _  r,  /3  +  °2 

+  4  +  4  V  2 


(6.6.3) 


For  large  positive  F  our  model  must  yield 


r  t 

As  S' r 


(6.6.4) 


and  this  requires  that 


As 


18 
3  +  a 


2  ’ 


(6.6.5) 


We  need  a  way  to  represent  a  for  an  arbitrary  orientation  of  the  coordinates.  The  structure  of  S“  is 
characterized  by 


(S')3 


(6.6.6) 


which  for  (6.6.2)  is 

1V-  30 -»*)/< 

[(3  +  aI)/2]3^2 


(6.6.7) 


W  ranges  from  ~l/\/6for  axisymmetric  expansion  to  l/\/6  for  axisymmetric  contraction.  Plane  strain  and 
shear  flow  correspond  to  IV  =0.  Using  (6.6.5)  to  express  a  in  terms  of  As,  and  then  in  turn  expressing  W 
in  terms  of  As,  we  find 


W  = 


(6.6.8) 


This  allows  us  to  determine  As  from  a  known  W.  The  relationship  between  them  is  shown  in  Fig.  6.6.1. 


0 


w 


1 

vB 

Figure  6.6.1  Variation  of  the  model  parameter  with  strain  type 


_1_ 


The  constant  .do  should  be  chosen  to  produce  the  proper  level  of  shear  stress  in  shear  flow,  for  that 
is  the  most  important  engineering  flow.  Shear  flow  can  be  represented  as  a  combination  of  rotation  and 
irrotational  strain.  Denoting  Vi>2  =  I'. 


fi.j  =  0  0 

\  0  0  0, 


(6.6.9a,  6) 


Hence,  for  shear  flow  IV  =  0,  d.e  -  3/\/2,  and  S  -  Y/\/2.  With  these  values,  do  =  23  produces  fc12  =  0.15 

at  Tr  =  12.7,  corresponding  to  the  homogeneous  shear  flow  experiments  of  Tavoularis  and  Oorrsin. 

We  now  have  an  anisotropy  model  that  is  always  realisable  for  all  types  of  strain,  and  has  the  right 
general  trend  of  btJ  with  5.1,  but  assumes  that  a  state  of  structural  equilibrium  has  been  attained.  In  order 
to  handle  transients,  we  propose  an  evolution  equation  for  btJ  that  would  give  (6.6.1)  as  its  equilibrium 
solution, 

K,  =  -Ci  [(i40  +  AsS’t)^,  +  S,)t]/t  (6.6.10) 

By  choosing  C i  =  4/15,  the  model  will  agree  with  the  initial  phase  of  rapid  distortion  of  isotropic  turbulence, 
and  the  rate  of  return  tj  isotropy  is  of  the  right  general  magnitude  for  linear  approximations.  Note  that  the 
model  correctly  predicts  no  change  in  the  ansiolropy  of  isotropic  turbulence  under  pure  rotation. 

For  many  engineering  problems  the  main  objective  of  the  turbulence  model  is  to  reveal  important  trends 
This  simple  anisotropy  model  would  make  the  important  stresses  change  in  the  right  general  way,  without 
becoming  unrealizable,  and  therefore  it  should  be  an  attractive  alternative  for  use  in  simple  two-equation 
turbulence  models.  Preliminary  studies  by  students  in  the  author’s  turbulence  class  support  this  conjecture. 


7.  NUMERICAL  SIMULATIONS  OF  TUBULENCE 
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7.1  Introduction 

Over  the  past  decade,  two  important  types  of  numerical  simulations  have  become  important.  The 
earlier  work  concentrated  on  large  eddy  simulations  (LES),  in  which  simple  models  are  used  for  the  small- 
scale  turbulence  and  a  realisation  of  the  large-scale  turbulence  is  computed.  The  underlying  idea  is  that  the 
structure  of  large  eddies  differ  greatly  fiom  flow  to  flow  (which  is  why  universal  models  are  elusive),  whereas 
the  small  eddies  are  more  universal  and  therefore  easier  to  model.  Large  eddy  simulations  have  provided 
important  information  for  turbulence  modeling,  and  there  is  now  great  interest  in  the  development  of  large 
eddy  simulations  as  a  tool  for  engineering  analysis.  A  prominent  program  in  this  direction  exists  in  FYance 
at  the  EDF. 

It  was  argued  that,  since  the  iatio  of  the  largest  to  the  smallest  scales  of  turbulence  varies  as  R '?J* 
(see  3.2.4),  it  would  never  be  practical  to  do  a  significant  simulation  of  all  the  important  turbulent  scales. 
However,  valid  direct  simulations  of  turbulent  flows  at  Rt  of  the  order  of  100-300  have  become  possible.  This 
is  the  range  of  turbulence  Reynolc  numbers  in  turbulent  shear  flows  with  Reynolds  numbers,  based  on  the 
layer  thickness  and  the  driving  mea  >  elocity  difference,  of  about  1000,  and  a  number  of  direct  simulations  of 
channel  flows  and  boundary  layer  flows  at  these  low  Reynolds  numbers  have  now  been  attained.  These  direct 
simulations  provide  an  important  new  tool  for  studying  turbulence,  particularly  because  they  yield  essentally 
any  data  that  one  might  desire.  Already  they  have  contributed  important  new  insight  into  turbulent  structure 
and  have  aided  advances  in  turbulence  modeling,  as  well  as  new  understanding  of  transition  physics. 

In  this  chapter  we  will  review  the  fundamentals  and  current  status  of  this  very  fast-moving  area  of 
research,  drawing  primarily  from  the  experience  of  the  large  group  working  in  this  area  at  the  NASA/Ames 
Research  Center  and  Stanford  University.  At  present  this  group  involves  about  ten  NASA  scientists,  three 
Stanford  Professors,  a  dozen  or  so  graduate  students,  and  some  post-doctoral  scholars  and  other  visitors, 
with  the  work  being  coordinated  by  the  joint  NASA/Stanford  Center  for  Turbulence  Research.  Some  of  the 
exciting  new  things  going  on  in  this  group  will  be  outlined,  with  details  being  left  for  the  authors  to  report 
for  themselves. 

7.2  Fundamentals  of  large  eddy  simulation 

In  LES  one  needs  a  way  to  define  the  large-scale  components  of  the  fields,  and  filtering  is  usually  used. 
The  filtered  field  f  is  defined  by 


/(x,t)  =  j  C(x,x';A)/(x',t)dV. 


(7.2.1) 


Here  G  is  a  filter  function,  which  determines  exactly  what  fraction  of  the  motion  is  defined  as  being  large 
scale,  and  A  is  a  filter  parameter  that  implements  this  choice.  The  filter  function  must  be  normalized  such 
that 

I  C(x,x';  A)d3x'  =  1  (7.2.2) 


for  all  x.  The  residual  field  /'  is  then  what  is  left  over  after  filtering, 

/(x,  t)  =  7(x,t)  +  /'(x,f). 
The  filtered  residual  field  is  not  zero  since 

7  *  7  T  *  o. 

Filtering  (7.2.3), 


(7.2.3) 


(7.2.4) 


/  =  /  +  /'  (7.2.5a) 

so  the  filtered  residual  field  can  be  expressed  in  terms  of  the  singly  and  doubly-filtered  resolved  fields, 


/'  =  /-/• 


(7.2.56) 


I  v; 


Thi«  proves  very  useful  in  modeling  the  residual  turbulence. 

In  homogeneous  turbulence  the  filter  must  be  of  the  form  <7(x  -  x';  A).  Then 

7(x,t)  =  j  G(x  -  x';  A)/(x',  t)d?x' 

has  the  Fourier  transform 

/(k.  I)  =  <5(*;A)/(k,0  (7.2.6) 

where  the  k  argument  of  G  is  the  magnitude  of  the  k  vector.  Several  filters  ahve 
cut-off  filter 

C  it  \k,  -  *'|  <  kc 
0  otherwise 

make  a  clean  separation  of  large  and  small  scales  in  spectral  space,  but  the  Gibbs 
Fourier  transform  make  the  physical-space  interpretation  undesirable.  Smoother 
with  the  Gaussian  filter, 

C(x-x';A)  =  M’  (7.2.8) 

where  A  is  a  constant  determined  by  the  normalization  and  depends  on  the  number  of  directions  in  which 
the  filter  is  applied.  The  Fourier  transform  of  the  Gaussian  filter  is  also  Gaussian, 

<J(k;A)  =  fle-k’A,/JV  (7.2.9) 

Filtering  is  more  of  a  problem  for  inhomogeneous  flows  .  The  most  satisfying  approach  is  to  use 
an  appropriate  set  of  expansion  functions  in  the  inhomogeneous  directions  and  then  to  define  the  filtered 
value  as  the  n-term  approximation.  However,  most  work  has  instead  used  finite-difference  methods  in  the 
inhomogeneous  directions  with  the  Gaussian  filter  in  the  homogeneous  directions,  and  taken  whatever  implicit 
filtering  is  provided  by  the  difference  scheme.  This  is  not  very  satisfying  because  it  leaves  the  computed  field 
ill  defined,  and  does  not  provide  a  systematic  way  for  estimation  of  the  energy  content  in  the  residual  field. 
This  is  one  of  the  unsatisfying  loose  ends  in  LES  that  needs  to  be  cleaned  up  by  some  good  research. 

The  evolution  equations  for  the  filtred  field  are  derived  by  filtering  the  Navier-Stoles  equations,  so  it  is 
important  that  the  filtering  definition  commute  with  differentiations  with  respect  to  both  time  and  space. 
The  Gaussian  filter  has  this  property,  and  so  homogeneous  turbulence  really  can  be  done  properly  with  LES 
using  the  Gaussian  filter.  If  p  -  p[t)  then  the  filtered  continuity  equation  is 


been  explored.  The  sharp 
(7.2.7) 

phenomena  in  the  inverse 
behavior  can  be  obtained 


p  +  pu,  =  0. 


(7.2.10) 


Subtracting  this  from  the  full  equation, 


0 


(7.2.11) 


so  the  residual  field  is  divergence-free,  and  if  p  =  constant  the  filtered  field  is  divergence-free.  Filtering 
the  momentum  equations,  assuming  p  is  constant  and  again  allowing  p  =  p(f),  the  equation  for  the  filtered 
velocity  field  is 


u,  +  (iquy),,-  =  --p„  -t-i/u,„y  . 
P 


(7.2.12) 


Representing  the  velocity  as  the  sum  of  filtered  and  residual  components, 


UjUy  =  u,uy  +  Rjj  (7.2.13) 

where  the  residual  stress  terms  are  _  __ 

Rij  =  u.Uy  +  u'u y  +  uju ' .  (7.2.14) 

In  LES  one  needs  to  model  A,y.  Given  this  model,  and  a  suitable  computer,  and  a  few  little  details  like 
boundary  and  initial  conditions,  tingle  realization t  of  turbulence  fields  can  be  generated.  In  homogeneous 
turbulence  this  appears  to  be  sufficient,  because  volume  averages  over  a  single  realization  seem  to  provide 
good  representations  for  ensemble  averages. 
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The  term  u,u,  does  not  need  to  be  modeled  because  it  can  be  computed  directly  by  filtering  the  product 
of  the  filtered  velocities.  This  is  easily  done  in  Fourier  space,  and  we  handle  this  term  this  way  now.  Our 
earlier  representation  of  this  in  terms  of  u^uy  +  L;y,  where  Lt]  was  the  Leonard  stress,  is  now  abandoned. 

It  may  be  noted  that  we  have  not  made  any  mention  of  numerical  methods  and  have  avoided  use  of 
the  term  sub-grid  scale  turbulence.  We  believe  that  it  is  important  to  cast  the  LES  equations  in  a  way  that 
is  independent  of  the  numerical  method,  and  would  lend  itself  to  purely  theoretical  analysis.  However,  in 
reality  the  filter  width  that  is  taken  is  related  to  the  computational  grid  employed.  The  results  depend  upon 
the  ratio  of  filter  width  to  mesh  width,  and  the  best  results  are  obtained  when  the  filter  width  is  twice  the 
mesh  width. 

7.S  Modeling  the  residual  stresses  in  large  eddy  simulation 

One  can  not  afford  a  very  complex  model  for  the  residual  stresses  in  LES.  Almost  all  of  the  work  to 
date  has  been  done  with  simple  algebraic  models,  although  there  have  been  some  explorations  with  simple 
one-equation  turbulence  models. 

It  is  useful  to  separate  into  isotropic  and  anisotropic  parts,  as  is  done  with  viscous  stresses, 

R>i  =  ^Rkkbij  +  T,j.  (7.3.1) 

The  isotropic  term  is  absorbed  with  the  filtered  pressure  by  writing 

P'  =  -p+'-Rkk  (7.3.2) 

p  3 

and  then  P“  replaces  p/p  and  TtJ  replaces  R,,  in  (7.2.12). 

An  important  element  of  most  LES  calculations  is  the  Smagorinski  model,  which  assumes  that  the 
residual  TJy  is  a  linear  function  of  the  anisotropic  strain  rate  imposed  by  the  filtered  field 

T.y  =  -2  vrSij  (7.3.3) 

where  up  is  an  eddy  viscosity  of  the  residual  field.  If  it  is  assumed  that  the  length  scale  of  the  dominant 
residual  eddies  is  the  filter  width,  and  that  the  time  scale  is  that  set  by  the  strain  rate  of  the  filtered  field, 
then 

i T  =  (Cs&)2\/smnSmn.  (7.3.4) 

The  coeffient  in  this  model  can  in  principle  be  evaluated  by  performing  direct  numerical  simulations  on  a  fine 
mesh  (say  1283),  then  filtering  this  data  to  a  coarse  mesh  (say  83)  to  define  the  filtered  and  residual  fields, 
and  then  comparing  the  model  with  the  residual  field  from  the  coarse  filtering.  Clark  et.  al.  (1979)  were 
the  first  to  emply  this  technique,  which  is  now  known  as  a  Clark  test.  For  isotropic  turbulence  the  results 
are  moderately  encouraging,  and  do  not  show  much  dependence  on  Reynolds  number,  a  value  of  about  0.12 
being  typical.  However,  when  this  test  is  applied  in  strained  and  sheared  flows,  essentially  no  correlation 
is  found  between  the  model  and  the  data.  The  model  simply  is  inadequate  under  these  more  interesting 
circumstances. 

An  important  advance  in  residual  stress  modeling  was  made  by  Bardina  (1985),  who  first  proposed  to 
model 

R,,  =  CD  (u,  ttj  -  Ujiiy) .  (7.3.5) 

The  basic  idea  was  to  characterise  the  stresses  of  the  residual  scales  as  being  similar  to  that  of  the  smallest 
resolvable  motions,  so  Bardina  called  this  the  scale  similarity  model.  By  itself  it  was  not  adequate  either, 
because  it  does  not  dissipate  sufficient  energy.  But  it  does  provide  energy  transfer  from  high  to  low  wavenum¬ 
bers,  and  effect  that  is  missing  in  the  Smagorinsky  model.  When  used  in  combination  with  the  Smagorinsky 
model  (the  Bardina  raized  model)  remarkably  good  results  are  obtained  in  the  Clark  tests,  with  the  same 
values  of  the  constant  yielding  correlations  between  predicted  and  actual  stresses  of  the  order  of  70%  for 
shear  flow,  irrotational  strain,  and  unstrained  flow!. 


The  value  of  the  constant  Cg  can  actually  be  deduced  from  a  simple  theoretical  argument.  If  one 
transform  to  new  coordinates  moving  linearly  with  respect  to  the  original  ones, 

i*  =  X{  -  Cjt  t'  =  t  u*  =  u,  +  Ci  (7.3.6a,  b) 


the  equations  of  motion  of  course  do  not  change  because  they  are  invariant  under  such  (Galelean)  trans¬ 
formations.  However,  individial  terms  in  the  equations  do  change  when  transformed.  For  the  filtering 
operation, 


U\  =  Ui  -f  Ci 


(7.3.7a,  6,  c) 


so  that  Rij  transforms  to 


Rij  =  u\u*'.  +  U*;.U*y  +  +C,U*J  +  Cyti-'  =  R’j  +  C,u'  +  Cyu'. 


(7.3.8) 


The  terms  modeling  JByy  should  transform  in  the  same  way;  the  Smagorinshy  model  is  invariant  under  the 
transformation,  and  hence  can  not  possibly  represent  all  of  Rij.  The  added  terms  of  the  Bardina  model 
(7.3.5)  transform  to 


Rij  =  Co  [u\u'y  -  U\«‘y  +Ci(u*y  -  u‘y)  +  Cy(u\  -  u\)].  (7.3.9) 

Using  (7.2.3)  and  (7.3.7b),  this  becomes 

Rij  =  R‘j  +  CB  (c,u^  +  cX) .  (7.3. 10) 

Comparing  (7.3.8)  and  (7.3.10),  it  is  evident  that  Cg  =  1.  Bardina  was  unaware  of  this  result  at  the  time 
he  did  his  numerical  work,  on  the  basis  of  which  he  recommended  a  value  of  1.05! 

In  recent  work  yet  to  be  published,  Piomelli  has  been  reexamining  LES  residual  modeling  using  the 
recent  direct  simulation  of  channel  flow  as  the  basis  for  Clark  tests,  also  carrying  out  LES  simulations  with 
various  models.  This  work  has  shed  some  new  light  on  LES  modeling,  which  can  be  summarized  as  follows. 
In  coarse  mesh  calculations  (say  163)  no  real  difference  is  observed  between  using  just  the  Smagorinsky 
model  and  the  Bardina  mixed  model,  and  the  results  in  general  reflect  the  coarseness  of  the  grid.  However, 
at  643  calculations  there  are  important  differences.  The  calculations  are  filtering  has  been  in  planes  parallel 
to  the  wall  only,  because  as  yet  we  do  not  really  have  any  good  way  to  do  explicit  filtering  in  directions  of 
inhomogeneity.  Piomelli  finds  that  the  choice  of  filter  function  is  important  in  determining  the  performance 
of  the  residual  turbulence  model.  The  filter  makes  its  appearance  in  the  calculations  when  the  term  u,uy  is 
calculated  by  filtering  the  product  of  the  computed  filtered  components.  If  the  Gaussian  filter  is  used  with 
the  Bardina  mixed  model,  very  good  results  are  obtained.  If  the  Gaussian  filter  is  used  with  the  Smagorinsky 
model,  very  poor  results  are  obtained.  But  if  the  Smagorinski  model  is  used  with  the  sharp  cut-off  filter,  fair 
results  are  obtained. 

The  inference  from  this  work  is  that  the  sharp  cut-off  filter  defines  a  clear  length  scale  for  the  residual 
turbulence,  whereas  the  Gaussian  filter  spreads  the  residual  scales  out  over  a  broader  range.  The  Bardina 
model  accounts  for  the  different  scales  in  the  residual  field  generated  by  the  Gaussian  filter.  On  the  other 
hand,  only  one  length  scale  is  carried  by  the  Smagorinsky  model,  and  therefore  this  model  can  not  account 
for  all  the  scales  filtered  by  the  Gaussian  filter. 

One  might  argue  that  the  turbulence  time  scale  in  the  Smagorinsky  viscosity  should  be  a  scale  appro¬ 
priate  to  the  residual  field.  In  isotropic  turbulence  the  strain  rate  of  the  resolved  field  sets  this  scale,  but  in 
inhomogeneous  flows  with  strong  mean  strain  rate  it  may  be  better  to  extract  the  time  scale  from  the  high 
wavenumber  end  of  the  resolved  field,  as  in  the  Bardina  model.  One  possible  approach  is  to  use  the  velocity 
scale  in  this  range,  _ 

ier  =  CA\J(ui,  -  uk)(uk  -  Ufc)- 
Another  approach  would  be  to  use  the  strain  rate, 


(7.3.11) 
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In  LES  one  probably  does  not  want  to  attempt  to  resolve  the  wall  region  of  boundary  layers,  and  so 
some  appropriate  wall  conditions  are  needed.  For  high  Reynolds  numbers,  it  is  through  this  condition  that 
the  viscosity  will  enter  the  problem.  The  main  thrust  of  Piomelli's  work  has  been  to  assess  various  proposals 
for  these  conditions.  At  this  writing  about  all  we  can  say  is  that  nothing  that  we  or  anyone  else  has  suggested 
shows  up  very  will  in  Clark  tests  against  the  direct  simulations  of  channel  flow.  However,  we  are  hopeful 
that  a  satisfactory  working  model  for  the  residual  wall  stress  will  be  found,  and  this  probably  will  draw 
upon  new  knowledge  about  the  structure  of  the  wall  region  that  is  currently  being  extracted  from  the  direct 
simulations. 


7.4  Insights  from  direct  simulations  of  homogeneous  turbulence 

Boundary  conditions  are  a  problem  in  turbulence  simulations.  The  problem  is  avoided  in  homogeneous 
turbulence  by  use  of  periodic  boundary  conditions.  The  resulting  turbulence  is  somewhat  artificial  in  that 
the  motion  on  opposite  sides  of  the  computational  domain  is  fully  correlated,  which  of  course  would  not  be 
the  case  in  a  real  turbulence  field.  One  must  select  a  computational  domain  large  enough  that  the  statistical 
correlations  at  separations  of  half  the  computational  domain  are  small,  and  when  this  is  done  the  statistical 
results  up  to  this  separation  seem  to  be  quite  like  those  of  real  turbulence. 

A  large  number  of  homogeneous  turbulence  simulations  have  been  carried  out  by  the  Ames/Stanford 
group,  almost  all  using  the  Rogallo  code.  This  program  uses  the  coordinate  transformation  (4.2.4),  and  as 
a  result  achieves  remarkable  robustness  in  runs  with  very  strong  deformation.  For  a  recent  description  of 
the  code  see  Lee  and  Reynolds  (1985).  Simulations  now  include  homogeneous  shear  flow  at  a  variety  of 
shear  rates,  many  cases  including  scalar  transport,  a  variety  of  irrotational  strain  flows,  return  to  isotropy 
following  various  strains,  some  rotation  cas^s.  Special  codes  have  handled  a  funny  type  of  homogeneous 
compressible  shear  flow  and  some  flow  compression  cases.  Meshes  ranging  from  643  to  2563  have  been  used, 
although  the  1283  cases  are  now  the  most  abundant. 

In  a  direct  simulation  one  must  capture  both  the  energy  at  large  scales  and  the  dissipation  at  small 
scales,  and  this  limits  the  calculations  to  relatively  low  Reynolds  numbers.  One  can  usually  tell  when  not 
enough  small-scales  have  been  captured  by  a  pile-up  of  energy  at  the  high  wavenumebr  end  of  the  spectrum. 
The  the  model  spectrum  (3.9.4)  can  be  used  to  estimate  the  fraction  of  energy  left  out  of  a  calculation  at 
any  given  Rt .  Typical  1283  calculations  miss  less  than  1%  of  the  turbulence  energy  at  Rt  =  50,  a  typical 
range  for  these  simulations. 

The  initial  turbulence  field  must  be  constructed  in  a  divergence-free  manner,  and  this  is  easily  done  with 
the  Fourier  representation.  The  spectrum  can  be  shaped  initially  and  scaled  to  contain  the  proper  energy 
for  a  target  Rt-  For  details  see  Lee  and  Reynolds  (1985). 


Figure  7.4.1  Spectra  for  relaxation  from  plane  strain 
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All  of  these  calculations  show  a  remarkable  amount  of  small-scale  anisotropy.  For  example,  Fig. 7. 4.1  shows 
one  of  Lee’s  spectra  during  relaxation  from  plane  strain,  with  the  different  lines  representing  different  com¬ 
ponents.  Note  that  the  anisotropy  persists  throughout  the  -S/3  range  of  the  spectrum.  We  investigated 
this  issue  of  small  scale  anisotropy  by  extending  the  measures  of  anisotropy  discussed  in  Chapter  6  to  the 
vorticity  and  dissipatin  fields.  The  vorticity  tensor  is  defined  as 


and  the  vortictly  anisotropy  tensor  is 


VH  = 


_  Vtj  -  u*8if/ 3 


The  dissipation  anisotropy  tensor  is  defined  by 

4/  =  D']  ~ n"6'’ 13 ■  (7-4.3) 

These  two  anisotropy  tensors  are  characterised  by  their  second  and  third  invariants,  defined  the  same  way 
as  those  for  the  Reynolds  stress  anisotropy  tensor  6,-y  (see  6.1.3).  Their  anisotropy  invariant  maps  are  the 
same  form  as  those  for  bl}  explained  in  section  6.1.  The  boundary  lines  are  the  same  for  the  6,y  and  d,j 
invariant  maps,  but  on  the  vorticity  invariant  map  the  two  axisymmetric  side  boundaries  are  reversed,  and 
the  uppermost  point  corresponding  to  one-dimensional  vorticity  corresponds  to  the  two-dimensional  velocity 
field. 

Fig.  (7.4.2)  shows  the  second  invariants  of  vorticity  and  velocity  during  relaxation  to  isotropy  from 
a  variety  of  different  strain  types.  The  trajectories  on  this  diagram  are  generally  down  and  then  to  the 
left.  Upon  the  removal  of  mean  strain  rate,  the  vorticity  anisotropy  relaxes  quickly  to  a  point,  and  then 
relaxes  slowly,  locked  on  to  the  anisotropy  of  the  Reynolds  stress!.  Moreover,  essentially  all  of  the  points 
showed  more  anisotropy  of  the  vorticity  thatn  of  the  Reynolds  stress!  These  are  astonishing  observations  to 
anyone  who  has  grown  up  with  the  idea  that  the  small  scales  become  isotropic  quickly,  compared  to  the  slow 
relaxation  of  the  scale  anisotropy. 

It  is  also  very  interesting  that  the  relationship  between  the  two  invariants  in  the  lock-on  phase  returning 
from  axisymmetric  expansion  is  quite  different  than  that  when  returning  from  axisymmetric  contraction. 
This  suggests  that  there  may  be  two  types  of  competing  structures  in  turbulence,  the  noodles  formed  by 
axisymmetric  contraction  and  the  pancakes  formed  by  axisymmetric  expansion,  and  that  perhaps  better 
turbulence  models  could  be  made  by  treating  these  structures  separately. 

We  have  mentioned  that  the  trajectories  for  return  to  isotropy  on  the  Reynolds  stress  invariant  map  are 
not  well  behaved,  which  casts  doubt  on  the  viability  of  modeling  the  slow  pressure  strain  and  dissipation 
anisotropy  terms  in  terms  of  t,y.  However,  those  on  the  vorticity  map  are  extremely  well  behaved.  Figure 
(7.4.3)  shows  these  trajectories,  which  are  well  fit  by  the  simple  model 


Vi,  =  -a  — Vi, 


where  a  depends  on  both  the  invariants  of  6,-y  and  o,y.  The  dissipation  anisotropy  trajectories  are  quite 
different,  but  they  too  are  very  well  behaved  and  can  be  modeled  quite  neatly.  For  details  see  Lee  and 
Reynolds  (1985). 

Upon  reflection,  the  requirement  that  the  vorticity  field  be  anisotropic  is  obvious  from  the  Biot-Savart 
law;  if  the  vorticity  spectrum  were  isotropic,  the  Reynolds  stress  spectrum  would  be  isotropic.  It  may  be 
that  explicit  consideration  of  this  anisotropy  in  turbulence  modeling  could  have  some  advantages.  We  have 
been  exploring  some  possibilities. 

In  another  recent  study,  Rogers  (1986)  has  examined  the  structure  of  homogeneous  turbulent  shear 
flow.  His  studies  reveal  that  hairpin  vortices  of  the  type  found  in  wall  boundary  layers  are  also  found  in 
homogeneous  turbulence.  However,  in  homogeneous  turbulence  there  are  both  “up"  and  “down"  hairpins, 
while  in  a  boundary  layer  one  sees  only  one  kind.  He  also  found  evidence  of  some  transverse  vortices  believed 
to  be  associated  with  the  weak  orientation  of  vorticity  caused  by  mean  rotation  (see  section  4.6). 
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Lee  has  extended  Rogers  work  to  (high)  shear  rates  and  Reynolds  numbers  comparable  with  the  viscous 
region  of  turbulent  boundary  layers.  Remarkably,  he  finds  long  longitudinal  streaks  that  familiar  objects  in 
the  wall  region,  with  transverse  spacings  that  scale  on  the  turbulent  stress  and  viscosity  in  exactly  the  same 
way  as  in  wall  boundary  layers.  This  work  suggests  that  it  is  the  high  shear  rate,  and  not  the  wall,  that 
produces  the  streaks!  This  would  be  good  news  for  modelers,  because  it  would  mean  that  models  based  on 
homogeneous  turbulence  might  have  far  more  to  do  with  boundary  layer  flows  than  one  might  think. 

Rogers  also  studied  scalar  transport  in  homogeneous  shear  flow  at  three  different  Prandtl  numbers 
There  are  three  interesting  situations  corresponding  to  an  (imposed)  means  scalar  gradient  in  each  direction. 
He  calculated  the  scalar  fields  for  all  three  cases  at  the  same  time  for  a  set  of  common  hydrodynamic 
simulations.  A  surprising  result,  actually  seen  in  experiments  by  Tavoularis  and  Corrsin,  is  that  some 
cross-gradient  scalar  fluxes  are  larger  than  the  flux  in  the  direction  of  the  mean  gradient!. 

Rogers  used  his  insight  about  the  hairpin  vortex  structures  and  the  transverse  vortices  to  explain  the 
mechanism  by  which  these  cross-gradient  transports  can  develop.  He  then  went  on  to  model  the  scalar 
flux  in  two  ways,  using  his  simulation  data  both  as  a  guide  in  the  modeling  and  as  the  basis  for  coefficient 
evaluation.  The  models  deal  with  an  anisotropic  diffusion  tensor  Z?,-y,  defined  by 

h,  =  W~6'  =  (7.4.5) 

where  S'  is  the  scalar  fluctuation  and  0O  is  a  mean  scalar  gradient.  The  diffusivity  tensor  could  be  calculated 
from  his  measurements,  and  is  found  to  be  inherently  non-symmetric.  However,  he  did  find  that  it  became 
antisymmetric  in  a  coordinate  system  that  is  aligned  with  the  principal  axes  of  the  Reynolds  stress.  This 
led  him  to  model  the  diffusion  tensor  in  the  form 

A j  —  Ci6,j  r  C'2  Aj  +  (7.4.6a, b) 


He  was  able  to  correlate  his  coefficients  with  Reynolds  and  Prandtl  numbers  to  within  about  20%. 

Rogers  made  another  model  assuming  that  the  scalar  flux  is  aligned  with  the  sum  of  the  mean  gradient 
terms  in  its  own  transport  equation,  and  thereby  obtained  a  model  of  comparable  accuracy  with  only  one 
free  coefficient.  This  model  is 


l  (  1.17\ 0  152  / 

-Ci)h,  +  hjV,„  +RijQ,j  =  0  Co  —  16.1  (  1  +  -jjj—  J  I  1  + 


(7.4.7a,  6) 


where  r  =  q2/e  and  Rt  =  94/(ice).  This  result  should  be  of  immediate  use  in  turbulence  modeling  for  both 
homogeneous  and  inhomogeneous  flows.  Rogers  has  recently  checked  this  model  against  direct  simulations 
of  turbulent  channel  flow  at  Pr  =  1  and  found  that  it  is  remarkably  accurate  for  the  flux  in  the  direction  of 
the  mean  temperature  gradient  and  within  about  20%  for  the  flux  perpendicular  to  the  mean  temperature 
gradient. 


7.5  Direct  simulations  of  spatially-developing  flows 


Some  of  the  most  exciting  work  at  present  are  the  boundary  layer  simulations  of  Spalart.  He  is  using 
a  clever  stretching  of  the  coordinate  system  that  enables  him  to  use  periodic  inflow-outflow  conditions  in  a 
growing  boundary  layer,  and  has  already  produced  results  about  the  structure  of  boundary  layers  in  pressure 
gradients  of  much  interest  to  experimentalists. 

In  order  to  simulate  more  general  turbulent  flows,  inflow  and  outflow  conditions  are  needed.  The  outflow 
problem  is  simpler  and  we  have  had  a  reasonable  solution  for  some  time.  The  inflow  problem  is  harder,  but 
we  have  recently  made  some  excellent  progress. 

Lowery  (1986)  simulated  the  spatially-developing  mixing  layer,  including  scalar  transport.  He  found 


that  a  soft  convective  outflow  condition, 


R  +  *K- 


(7.5.1) 


applied  to  the  velocity  components  and  scalar  worked  quite  well,  with  minumum  upstream  influence.  The 
convection  velocity  Ue  was  taken  as  the  average  of  the  two  free  stream  speeds.  At  the  inlet  he  forced  the 
flow  with  a  combination  of  fundamental  and  two  subharmonics  of  a  dominant  instability  of  the  inlet  layer 


(tanh  profile),  because  the  layer  was  forced,  it  responded  like  a  forced  layer,  with  parings  occuring  cyclicly 
at  froien  locations.  And,  the  layer  grows  not  linearly,  as  do  natural  layers,  but  by  leaps  and  bounds,  as  do 
forced  mixing  layers  in  the  laboratory. 

It  has  been  asked  if  the  mixing  layer  is  absolutely  unstable,  in  which  case  if  the  forcing  is  stopped  after 
large  disturbances  have  developed  downstream  the  layer  should  continue  to  remain.  When  Lowery  terminated 
forcing,  the  initial  region  of  the  layer  relaminarised,  suggesting  that  the  instability  was  convective,  but  midway 
down  the  flow  the  turbulence  never  went  away,  and  by  the  exit  the  flow  was  quite  turbulent.  His  calculation 
did  not  include  the  splitter  plate,  which  undoubtedly  plays  a  role  in  promoting  absolute  instability,  so  the 
matter  is  not  really  resolved.  Lowery  also  studied  the  growth  of  three-dimensional  disturbances  in  the  layer, 
finding  that  they  grew  to  scales  and  structures  characteristic  of  the  braid  region  of  the  mixing  layer. 

Ongoing  extensions  of  our  mixing  layer  simulation  work  by  Sandham  involve  the  use  of  random  jitter 
of  the  forcing  to  simulate  more  natural  turbulent  inflow  condition.  This  produces  the  linear  growth  seen 
in  natural  experimental  layers,  at  growth  rates  in  excellent  agreement  with  experiemnts.  The  resulting 
statistical  quantisers,  including  the  scalar  pdf,  are  much  more  like  those  measured  for  natural  layers.  It 
now  seems  that  this  will  be  quite  an  acceptable  method  for  generating  relativelty  simple  yet  effective  inflow 
conditons  for  direct  numerical  simulations  of  turbulence. 

Current  work  is  concentrating  on  extensions  to  compressible  mixing  layers,  the  goal  being  to  use  these 
direct  simulations  as  the  basis  for  building  better  turbulence  models  for  supersonic  flows,  including  combus¬ 
tion,  both  for  use  in  LES  and  in  simpler  turbulence  models.  There  is  a  growing  group  at  Ames,  involving 
Rogers,  Moser  and  others,  beginning  to  work  very  seriously  on  turbulent  combustion  simulations.  It  seems 
safe  to  forecast  that  a  decade  from  now  the  capabilities  for  know  much  more  about  the  modeling  and  simula¬ 
tion  of  these  and  flows  of  technical  interest  will  be  considerably  advanced,  and  students  who  have  mastered 
these  notes  should  be  ready  to  begin  the  exciting  work  ahead  in  this  area. 
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